Population analysis

In this script, we load the population data and compute summary statistics and visualisations.

Clear workspace and add repositories to path

clear all
% path_repo = fileparts(mfilename('fullpath'));
path_repo = "<YOUR_LOCAL_PATH>/Thiele-attention-gratc-LIP-pipette/repositories/attention-parietal-dopamine/analyses/";
% https://github.com/jochemvankempen/attention-parietal dopamine
addpath(genpath(fullfile(path_repo,'/../../attention-parietal-dopamine')))
% https://github.com/jochemvankempen/plotj
addpath(genpath(fullfile(path_repo,'/../../plotj')))
% https://github.com/jochemvankempen/gain-variability
addpath(genpath(fullfile(path_repo,'/../../gain-variability')))

Directories

Set data directories

subjects = {'W','J'};
% data path, relative to repository
path_data = fullfile(path_repo,'/../../../data/processed');
path_population = regexprep(path_data,'processed','population');
if ~isfolder(path_population)
    mkdir(path_population)
end
fprintf('Load data from: %s\n', path_data)
Load data from: /Users/jochemvankempen/NCL/gratc_DA/Thiele-attention-gratc-LIP-pipette/repositories/attention-parietal-dopamine/analyses/../../../data/processed
fprintf('Save data to: %s\n', path_population)
Save data to: /Users/jochemvankempen/NCL/gratc_DA/Thiele-attention-gratc-LIP-pipette/repositories/attention-parietal-dopamine/analyses/../../../data/population

Param

set some parameters that we'll use

label_drug = {'Dopamine','SCH23390'};
label_drug_ext = {'',' (D1-antagonist)'};
label_attention = {'Attend RF','Attend away'};
label_drug_onoff = {'Drug on','Drug off'};
label_unitclass = {'Narrow','Broad'};
markersize = 15; % scatterplot markersize
plot_conventions = 'Nature'; % follow figure conventions from this journal
% condition indices
idx_attention_cond = [1 0 0 1 0 0]; % cond_num mapping to attend RF/away in trialdata
idx_attention_roc = [1 1 0 1 1 0]; % see condition comparisons in spike_rate_ROC.m
% select units selective for task feature, e.g. attention or drug selective
% selectivity_criterium = 'none';
% selectivity_criterium = 'dru';
% selectivity_criterium = 'dru';
% selectivity_criterium = 'att';
selectivity_criterium = 'att&dru';

Get recordinglist

Compile list of recordings from recinfo.mat in each session folder

recordinglist = get_recordinglist(subjects, path_data);
fprintf('Loaded %d recordings \n', height(recordinglist))
Loaded 47 recordings 

Get population data

Compile population data matrices/tables

[rate_summary, unitlist] = get_population_data(recordinglist, 'spike_rate_summary', path_data, 'dim', 'block_change(1-Inf)');
rate_PSTH = get_population_data(recordinglist, 'spike_rate_PSTH', path_data, 'dim');
rate_ANOVA = get_population_data(recordinglist, 'spike_rate_ANOVA', path_data, 'dim');
rate_ROC = get_population_data(recordinglist, 'spike_rate_ROC', path_data, 'dim', 'block_change(1-Inf)');
waveform = get_population_data(recordinglist, 'waveform', path_data);
RT = get_population_data(recordinglist, 'RT_drug_modulation', path_data, 'stim');

Extract param

num_rec = size(RT.RT,1);
[num_unit, num_drug, num_cond] = size(rate_summary.rate);
fprintf('Loaded %d units \n',num_unit)
Loaded 88 units 

Classify cells as narrow or broad spiking

waveform_cutoff = 250;
unit_class = cell(num_unit,1);
unit_class(waveform.peak_to_trough_time <= waveform_cutoff) = {'Narrow'};
unit_class(waveform.peak_to_trough_time > waveform_cutoff) = {'Broad'};
% fill and define as categorical variable
unitlist.peak_to_trough_time = waveform.peak_to_trough_time;
unitlist.unit_class = categorical(unit_class, 'Ordinal', false);
unitlist.unit_class = reordercats(unitlist.unit_class, label_unitclass);
unitlist.waveform = waveform.waveform;
waveform_time = waveform.time;
clear waveform unit_class

Append unitlist

% rename units in sequential order to be able to join to later tables
unitlist.unit = categorical((1:num_unit)', 'Ordinal', false);
% average over the ejection currents from both barrels
unitlist.EjectCurrent = nanmean(unitlist.EjectCurrent,2);
unitlist.EjectCurrent_centered = unitlist.EjectCurrent - mean(unique(unitlist.EjectCurrent));
% append/join with selectivity/ANOVA table
unitlist = join(unitlist, rate_ANOVA.selectivity, 'Keys', 'unit');
% define categorical variables
unitlist.Drug = categorical(cellstr(unitlist.Drug), 'Ordinal', false);
unitlist.Subject = categorical(cellstr(unitlist.Subject), 'Ordinal', false);
% define datetime variables
unitlist.Date = datetime(unitlist.Date, 'Format', "yyyy-MM-dd");

Create population tables: attention

table for attention conditions, including drug on/off

% number of elements in table
t_numel = num_unit*num_drug*2;
% compute attention averages
t_rate(:,:,1) = mean(rate_summary.rate(:,:,idx_attention_cond==1),3);
t_rate(:,:,2) = mean(rate_summary.rate(:,:,idx_attention_cond==0),3);
t_FF(:,:,1) = mean(rate_summary.FF(:,:,idx_attention_cond==1),3);
t_FF(:,:,2) = mean(rate_summary.FF(:,:,idx_attention_cond==0),3);
t_PSTH(:,:,1,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==1,:),3);
t_PSTH(:,:,2,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==0,:),3);
t_gain = arrayfun(@(x)(x.rhat),rate_ROC.gv,'UniformOutput',true);
t_gain_log = log(t_gain);
t_gv = rate_ROC.gv;
% use idx_attention_roc for ROC and MI, not that this is different from the
% other attention index.
t_roc_attend(:,:,1) = mean(rate_ROC.roc_attend(:,:,idx_attention_roc==1),3);
t_roc_attend(:,:,2) = mean(rate_ROC.roc_attend(:,:,idx_attention_roc==0),3);
t_mi_attend(:,:,1) = mean(rate_ROC.mi_attend(:,:,idx_attention_roc==1),3);
t_mi_attend(:,:,2) = mean(rate_ROC.mi_attend(:,:,idx_attention_roc==0),3);
% create metadata columns
[t_att] = deal(zeros(size(t_rate)));
t_att(:,:,1) = 1;
t_unit = categorical((1:num_unit)', 'Ordinal', false);
t_unit = repmat(t_unit, [1, num_drug, 2]);
t_drug_on = [0 1];
t_drug_on = repmat(t_drug_on, [num_unit, 1, 2]);
t_unit      = reshape(t_unit, [t_numel, 1]);
t_drug_on   = reshape(t_drug_on, [t_numel, 1]);
t_att       = reshape(t_att, [t_numel, 1]);
[t_att_label,t_drug_on_label] = deal(cell(length(t_att),1));
t_att_label(t_att==1) = {'Attend RF'};
t_att_label(t_att==0) = {'Attend away'};
t_drug_on_label(t_drug_on==0) = {'Drug off'};
t_drug_on_label(t_drug_on==1) = {'Drug on'};
% reshape matrices to single columns
t_rate          = reshape(t_rate, [t_numel, 1]);
t_FF            = reshape(t_FF, [t_numel, 1]);
t_roc_attend    = reshape(t_roc_attend, [t_numel, 1]);
t_mi_attend     = reshape(t_mi_attend, [t_numel, 1]);
t_gain          = reshape(t_gain, [t_numel, 1]);
t_gain_log      = reshape(t_gain_log, [t_numel, 1]);
t_gv            = reshape(t_gv, [t_numel, 1]);
t_PSTH          = reshape(t_PSTH, [t_numel, length(rate_PSTH.time)]);
% build table
att_table = table(...
    t_unit, ...
    t_drug_on_label, t_att_label, ...
    t_rate, t_FF, ... % rate, FF
    t_roc_attend, t_mi_attend, t_gain, t_gain_log, t_gv, ... % roc, MI, gain, gain_log, GainVariability
    t_PSTH, ... % PSTH
    'VariableNames', {'unit', 'drug_on', 'attention', 'rate', 'FF', 'roc', 'MI', 'gain', 'gain_log', 'gv', 'PSTH'});
% join tables
att_table = join(att_table, unitlist, 'Keys', 'unit');
% define categorical variables
att_table.drug_on = categorical(att_table.drug_on, 'Ordinal', false);
att_table.attention = categorical(att_table.attention, 'Ordinal', false);
% reorder categories to label_* order
att_table.drug_on = reordercats(att_table.drug_on, label_drug_onoff);
att_table.attention = reordercats(att_table.attention, label_attention);
% define combined classes
att_table.attention_drug_on = att_table.attention .* att_table.drug_on;
att_table.attention_drug_on = categorical(regexprep(cellstr(att_table.attention_drug_on), ' Drug', ', Drug')); % add a comma in between classes
% clear temporary variables
clear t_*
% print head
fprintf('att_table:\n')
att_table:
disp(head(att_table))
    unit    drug_on     attention     rate       FF        roc         MI         gain      gain_log             gv                   PSTH          Subject    Task        Date         Drug      EjectCurrent    Weight    peak_to_trough_time    unit_class       waveform       EjectCurrent_centered      s_stim        s_att         s_dru         s_dir       s_att*dru    s_att*dir    s_dru*dir    s_att*dru*dir     attention_drug_on 
    ____    ________    _________    ______    ______    _______    ________    ________    ________    _____________________    _______________    _______    _____    __________    ________    ____________    ______    ___________________    __________    ______________    _____________________    __________    __________    __________    __________    _________    _________    _________    _____________    ___________________

     1      Drug off    Attend RF    8.2778    2.1547     0.6376     0.12527     0.14343    -1.9419     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-05-21    Dopamine         25          9.29            345.6             Broad       [1×187 double]           -34.25                0.6978     0.0023999       0.57272       0.83275     0.55229      0.52515      0.30834        0.64143       Attend RF, Drug off
     2      Drug off    Attend RF    8.8506    1.4285    0.70263     0.17977         NaN        NaN     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-05-21    Dopamine         25          9.29            205.2             Narrow      [1×187 double]           -34.25               0.79882    5.7922e-14       0.16206       0.25224    0.027989      0.54777      0.21153        0.20304       Attend RF, Drug off
     3      Drug off    Attend RF    23.394    3.7531    0.66866     0.12168     0.11193    -2.1899     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-05-21    Dopamine         25          9.29            388.8             Broad       [1×187 double]           -34.25               0.05165    4.2769e-08       0.45615       0.31594    0.054758      0.42777      0.44512        0.26137       Attend RF, Drug off
     4      Drug off    Attend RF    11.453    1.8759    0.64466     0.11753    0.065949    -2.7189     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-05-28    Dopamine         60           NaN            394.2             Broad       [1×187 double]             0.75               0.29917    1.5365e-08       0.98854       0.90503      0.1141      0.97673      0.61297        0.38682       Attend RF, Drug off
     5      Drug off    Attend RF    28.963    5.2381    0.61176    0.088129     0.18975     -1.662     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-05-29    Dopamine         80           NaN            183.6             Narrow      [1×187 double]            20.75            2.0498e-22    6.6763e-07     0.0087419       0.88981     0.43338      0.71566      0.67429        0.44351       Attend RF, Drug off
     6      Drug off    Attend RF    44.259    7.4757    0.73118     0.19062     0.13876     -1.975     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-06-03    Dopamine         75             9            302.4             Broad       [1×187 double]            15.75            5.3262e-45    2.8411e-14    0.00011895       0.74406     0.24956      0.68049       0.2666         0.1016       Attend RF, Drug off
     7      Drug off    Attend RF    24.173    2.6038    0.73431      0.1636    0.065311    -2.7286     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-06-03    Dopamine         75             9            329.4             Broad       [1×187 double]            15.75            1.8937e-29    2.2529e-13      0.011994    0.00099457     0.80152     0.022988      0.93706        0.83731       Attend RF, Drug off
     8      Drug off    Attend RF    12.152    5.0028    0.57967     0.14417     0.33173    -1.1034     [1×1 GainVariability]    [1×1251 double]       W       gratc    2015-06-04    Dopamine         90             8            291.6             Broad       [1×187 double]            30.75            7.4239e-14    0.00078848      0.075078      0.043667     0.91433      0.34239      0.58849        0.20115       Attend RF, Drug off
% save to csv file for R analysis
table_selectivity = get_unit_selectivity(att_table, selectivity_criterium);
R_table = att_table(table_selectivity,:);
R_table = removevars(R_table,{'PSTH','waveform'});
writetable(R_table, fullfile(path_population,sprintf('population_attention_%s.csv',selectivity_criterium)))

drug table

table for analysis results that are computed for each attention condition, across drug conditions. There are thus no separate drug_on conditions in this table.

% number of elements in table
t_numel = num_unit*2;
% compute attention averages
t_roc_drug(:,1) = mean(rate_ROC.roc_drug(:,idx_attention_cond==1),2);
t_roc_drug(:,2) = mean(rate_ROC.roc_drug(:,idx_attention_cond==0),2);
t_mi_drug(:,1) = mean(rate_ROC.mi_drug(:,idx_attention_cond==1),2);
t_mi_drug(:,2) = mean(rate_ROC.mi_drug(:,idx_attention_cond==0),2);
t_roc_drug = reshape(t_roc_drug, [t_numel, 1]);
t_mi_drug = reshape(t_mi_drug, [t_numel, 1]);
% create metadata columns
t_unit = categorical((1:num_unit)', 'Ordinal', false);
t_unit = repmat(t_unit, [2, 1]);
t_cond = repmat(1:2,[num_unit 1]);
t_cond_label = cell(size(t_cond));
t_cond_label(:,1) = {'Attend RF'};
t_cond_label(:,2) = {'Attend away'};
t_cond_label = reshape(t_cond_label, [t_numel, 1]);
% build table
drug_table = table(...
    t_unit, ...
    t_cond_label, ...
    t_roc_drug, t_mi_drug, ... % roc, MI
    'VariableNames', {'unit', 'attention', 'roc', 'MI'});
% join tables
drug_table = join(drug_table, unitlist, 'Keys', 'unit');
% define categorical variables
drug_table.attention = categorical(drug_table.attention, 'Ordinal', false);
% reorder categories to label_* order
drug_table.attention = reordercats(drug_table.attention, label_attention);
% clear temporary variables
clear t_* rate_*
% print head
fprintf('drug_table:\n')
drug_table:
disp(head(drug_table))
    unit    attention      roc          MI        Subject    Task        Date         Drug      EjectCurrent    Weight    peak_to_trough_time    unit_class       waveform       EjectCurrent_centered      s_stim        s_att         s_dru         s_dir       s_att*dru    s_att*dir    s_dru*dir    s_att*dru*dir
    ____    _________    _______    __________    _______    _____    __________    ________    ____________    ______    ___________________    __________    ______________    _____________________    __________    __________    __________    __________    _________    _________    _________    _____________

     1      Attend RF    0.48333    0.00010785       W       gratc    2015-05-21    Dopamine         25          9.29            345.6             Broad       [1×187 double]           -34.25                0.6978     0.0023999       0.57272       0.83275     0.55229      0.52515      0.30834        0.64143   
     2      Attend RF    0.60052      0.085384       W       gratc    2015-05-21    Dopamine         25          9.29            205.2             Narrow      [1×187 double]           -34.25               0.79882    5.7922e-14       0.16206       0.25224    0.027989      0.54777      0.21153        0.20304   
     3      Attend RF    0.60693      0.068293       W       gratc    2015-05-21    Dopamine         25          9.29            388.8             Broad       [1×187 double]           -34.25               0.05165    4.2769e-08       0.45615       0.31594    0.054758      0.42777      0.44512        0.26137   
     4      Attend RF    0.57284      0.038441       W       gratc    2015-05-28    Dopamine         60           NaN            394.2             Broad       [1×187 double]             0.75               0.29917    1.5365e-08       0.98854       0.90503      0.1141      0.97673      0.61297        0.38682   
     5      Attend RF    0.44499     -0.034453       W       gratc    2015-05-29    Dopamine         80           NaN            183.6             Narrow      [1×187 double]            20.75            2.0498e-22    6.6763e-07     0.0087419       0.88981     0.43338      0.71566      0.67429        0.44351   
     6      Attend RF    0.36805     -0.097646       W       gratc    2015-06-03    Dopamine         75             9            302.4             Broad       [1×187 double]            15.75            5.3262e-45    2.8411e-14    0.00011895       0.74406     0.24956      0.68049       0.2666         0.1016   
     7      Attend RF    0.39486     -0.053193       W       gratc    2015-06-03    Dopamine         75             9            329.4             Broad       [1×187 double]            15.75            1.8937e-29    2.2529e-13      0.011994    0.00099457     0.80152     0.022988      0.93706        0.83731   
     8      Attend RF    0.48523      -0.08402       W       gratc    2015-06-04    Dopamine         90             8            291.6             Broad       [1×187 double]            30.75            7.4239e-14    0.00078848      0.075078      0.043667     0.91433      0.34239      0.58849        0.20115   

Print summary of data

Here we list the proportion of units that were selective for various

% all units
fprintf('UNITS\n')
UNITS
fprintf('\tRecorded %d units\n', height(unitlist))
	Recorded 88 units
% units per subject
fprintf('SUBJECTS\n')
SUBJECTS
M1 = unitlist.Subject==subjects{1};
M2 = unitlist.Subject==subjects{2};
fprintf('\tRecorded %d/%d units from monkey 1/2\n', length(find(M1)), length(find(M2)))
	Recorded 64/24 units from monkey 1/2
% unit selectivity (visual/attention)
fprintf('SELECTIVITY\n')
SELECTIVITY
unit_selectivity = get_unit_selectivity(unitlist, 'visual');
selective_visual = length(find(unit_selectivity));
fprintf('\t%d units with visual response (%1.1f%%)\n', selective_visual, selective_visual/num_unit*100)
	76 units with visual response (86.4%)
unit_selectivity = get_unit_selectivity(unitlist, 'att');
selective_attention = length(find(unit_selectivity));
fprintf('\t%d units with attentional modulation (%1.1f%%)\n', selective_attention, selective_attention/num_unit*100)
	74 units with attentional modulation (84.1%)
% units and unit selectivity per drug
fprintf('DRUG\n')
DRUG
selective_drug = 0;
selective_attdrug = 0;
selective_interaction = 0;
for idrug = 1:length(label_drug)
    
    num_drug_units = unitlist.Drug==label_drug(idrug); % total units recorded with this drug
    
    unit_selectivity = get_unit_selectivity(unitlist, 'dru',{'drug',label_drug(idrug)});
    
    fprintf('\tRecorded %d/%d (%d) units from monkey (1/2) with %s\n', length(find(M1 & num_drug_units)), length(find(M2 & num_drug_units)), length(find(num_drug_units)), label_drug{idrug})
    
    selective_drug = selective_drug+length(find(unit_selectivity));
    fprintf('\t%d units modulated by %s (%1.1f%%)\n', length(find(unit_selectivity)), label_drug{idrug}, length(find(unit_selectivity))/length(find(num_drug_units))*100)
    
    unit_selectivity = get_unit_selectivity(unitlist, 'att&dru',{'drug',label_drug(idrug)});
    selective_attdrug = selective_attdrug+length(find(unit_selectivity));
    fprintf('\t%d units modulated by attention and drug %s (%1.1f%%)\n', length(find(unit_selectivity)), label_drug{idrug}, length(find(unit_selectivity))/length(find(num_drug_units))*100)
    
    unit_selectivity = get_unit_selectivity(unitlist, 'att*dru',{'drug',label_drug(idrug)});
    selective_interaction = selective_interaction+length(find(unit_selectivity));
    fprintf('\t%d units interaction with attention %s (%1.1f%%)\n', length(find(unit_selectivity)), label_drug{idrug}, length(find(unit_selectivity))/length(find(num_drug_units))*100)
end
	Recorded 41/18 (59) units from monkey (1/2) with Dopamine
	36 units modulated by Dopamine (61.0%)
	31 units modulated by attention and drug Dopamine (52.5%)
	4 units interaction with attention Dopamine (6.8%)
	Recorded 23/6 (29) units from monkey (1/2) with SCH23390
	14 units modulated by SCH23390 (48.3%)
	14 units modulated by attention and drug SCH23390 (48.3%)
	5 units interaction with attention SCH23390 (17.2%)
% unit summary
fprintf('SUMMARY\n')
SUMMARY
fprintf('\t%d total units modulated by attention  (%1.1f%%)\n', selective_attention, selective_attention/num_unit*100)
	74 total units modulated by attention  (84.1%)
fprintf('\t%d total units modulated by drug application  (%1.1f%%)\n', selective_drug, selective_drug/num_unit*100)
	50 total units modulated by drug application  (56.8%)
fprintf('\t%d total units modulated by drug application and attention  (%1.1f%%)\n', selective_attdrug, selective_attdrug/num_unit*100)
	45 total units modulated by drug application and attention  (51.1%)
fprintf('\t%d total units interaction attention by drug application  (%1.1f%%)\n', selective_interaction, selective_interaction/num_unit*100)
	9 total units interaction attention by drug application  (10.2%)
%

Plot population histogram and example unit

[colors, labels] = get_colors('att_drug');

colors = reshape(colors, [4,3]);

labels = reshape(labels, [4,1]);

linestyle = {'-','-','-','-'};

events = {'stim', 'cue', 'dim'};

label_event = {'Stimulus onset', 'Cue onset', 'First-dimming'};

xlabel_event = {{'Time from','stimulus onset (ms)'}, {'Time from','cue onset (ms)'}, {'Time from','first-dimming (ms)'}};

xlim2use = {[-200 750], [-100 750], [-1000 200]};

ylim2use = [25 90];

ncol = length(label_event) + 2;

nrow = length(label_drug);

idx_subplot = {1,2,3,6,7,8,[4 5 9 10]};

[fH, fSet] = plotj_initFig('width', 26, 'height', 11, 'Journal',plot_conventions);

fSet.subplotGap2 = fSet.subplotGap .* [1.3 0.30];

iplot = 0;

for irow = 1:nrow

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(irow)});

clear h

for ievent = 1:length(label_event)

iplot = iplot+1;

rate_PSTH = get_population_data(recordinglist, 'spike_rate_PSTH', path_data, events{ievent});

times = rate_PSTH.time;

% convert samples to att_table format

t_numel = num_unit*num_drug*2;

clear t_PSTH

t_PSTH(:,:,1,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==1,:),3);

t_PSTH(:,:,2,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==0,:),3);

t_PSTH = reshape(t_PSTH, [t_numel, length(times)]);

% initialise axes

subtightplot(nrow, ncol, idx_subplot{iplot}, fSet.subplotGap2, fSet.subplotMargin, fSet.subplotMargin)

if ievent==1 && irow==1

plotj_initAx(fSet, 'axlabel', irow, 'axlabelDisplacement', [0.07, 0.00]);

else

plotj_initAx(fSet);

end

if ievent==(ceil(length(label_event)/2))

ht = title([label_drug{irow} label_drug_ext{irow}], 'FontSize', fSet.Fontsize_title);

% ht.Position = ht.Position + [0.05 0 0];

end

if (ievent~=1)

axh = get(gca);

axh.YTick = [];

axh.YAxis.Visible = 'off'; % remove y-axis

% axes('Color','none','YColor','none')

% set(gca,'xtick',[])

end

hold on

% plot each attention and drug condition

icol = 0;

conditions = categories(att_table.attention_drug_on);

for iatt_drug = 1:length(conditions)

idx_att_drug = att_table.attention_drug_on==labels{iatt_drug};

idx = idx_unit & idx_att_drug;

icol = icol+1;

h(icol) = boundedline(times, ...

squeeze(mean(t_PSTH(idx,:),1)), ...

squeeze(std(t_PSTH(idx,:),[],1))/sqrt(length(find(idx))), ...

'cmap', colors(iatt_drug,:), 'alpha');

set(h(icol), 'LineWidth', fSet.LineWidth)

set(h(icol), 'LineStyle', linestyle{iatt_drug})

end

xlim(xlim2use{ievent})

ylim(ylim2use)

if ievent==1

ylabel('Normalized firing rate (%)', 'FontSize', fSet.Fontsize_text)

text(200, 30, sprintf('n=%d', length(find(idx))))

elseif ievent==3

if 0%irow==1

hleg = legend(h, labels, 'FontSize', fSet.Fontsize_text, 'Box', 'Off');

hleg.Position = hleg.Position + [0.05 0.15 0 0];

end

if irow==nrow

xlabel(xlabel_event{ievent}, 'FontSize', fSet.Fontsize_text)

end

hleg = legend(h, labels, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'Location', 'NorthWest');

hleg.Position = hleg.Position + [-0.04 0.53 0 0];

% plot example cell, plot raster for example cell, only drug/attention

% plot param

iplot = iplot+1;

axH = subtightplot(nrow, ncol, idx_subplot{iplot}, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

hold on

plotj_initAx(fSet, 'axlabel', 2, 'axlabelDisplacement', [0.07, 0.00]);

yoffset = 100;

[colors] = get_colors('att_drug');

% load data for single unit

event2plot = 3;

rec2plot = {'W','2015-10-01'};

unit2plot = [41 1];

idx_unit = unitlist.Subject==rec2plot{1} ...

& unitlist.Date==rec2plot{2} ...

& unitlist.unit==categorical(unit2plot(1));

loadfilename = fullfile(path_data, string(unitlist.Subject(idx_unit)), string(unitlist.Date(idx_unit)), 'unit.mat');

unit = load(loadfilename);

loadfilename = fullfile(path_data, string(unitlist.Subject(idx_unit)), string(unitlist.Date(idx_unit)), 'trialdata.mat');

load(loadfilename);

% [trialdata, unit] = remove_excluded_trials(trialdata, unit, unit2plot);

% select only trial window for which this unit has spikes

sps = unit.Dim1Align;

% raster for each condition

idx_attend = idx_attention_cond([trialdata.cond_num])';

idx_drug = [trialdata.drug]' + 1;

ymax = 0;

icond = 0;

for iatt = 1:length(unique(idx_attention_cond))

for idrug = 1:2

icond = icond+1;

% get trial indices

trial_index = (...

idx_drug==idrug) ...

& abs(idx_attend-2)==iatt;

% compute raster

[xspikes, yspikes] = spike_raster( sps(:,trial_index), xlim2use{event2plot}, unit2plot(2));

% get PSTH

x = rate_PSTH.time;

y = squeeze(mean(rate_PSTH.samples(idx_unit,idrug,abs(idx_attention_cond-2)==iatt,:),3));

% plot

h(icol) = plot(x, y, 'Color', colors(iatt,idrug,:), 'linew', fSet.LineWidth, 'LineStyle', linestyle{icond});

plot(xspikes, yspikes+yoffset+ymax, 'Color', colors(iatt,idrug,:), 'linew', 1.5)

% set y-offset for next iteration

ymax = ymax + max(yspikes(:));

end

plot([-500 0], [10 10], 'Color', [0.5 0.5 0.5], 'linew', 5)

set(gca,'YTick',[0 yoffset],'YTicklabel',[0 round(rate_PSTH.maxhist(idx_unit)*1000)]);

xlabel('Time from first-dimming (ms)', 'FontSize', fSet.Fontsize_text)

ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)

xlim([-800 200])

ylim2use = yoffset+ymax+80;

ylim([0 ylim2use])

% set text in plot

p = [unitlist.s_att(idx_unit) unitlist.s_dru(idx_unit) unitlist.("s_att*dru")(idx_unit)];

[pstring,starstring] = get_significance_strings(p, 'rounded', 0, 'factorstring', {'attention', 'drug', 'interaction'});

text(axH, -750, ylim2use*0.95, pstring, 'FontSize', fSet.Fontsize_text)

Violin population plots, mean-var relationship, auroc

data2plot = {'rate','FF','gain_log'};

colors = get_colors('spikewidth');

markerstyle = {'o','o'};

[colors_violin, labels_violin] = get_colors('att_drug');

colors_violin = reshape(colors_violin, [4,3]);

labels_violin = reshape(labels_violin, [4,1]);

% plotting param

ncol = length(data2plot) + 1 + 1; %

nrow = 2; % drugs

iplot = 0;

plotlabel_offset = 0;

xval_gain = 1:500;

[fH, fSet] = plotj_initFig('figNum', 2, 'width', 26, 'height', 13, 'Journal',plot_conventions);

fSet.subplotGap = [0.1, 0.06];

fSet.subplotMargin = fSet.subplotMargin / 2;

for idrug = 1:num_drug

fprintf('Testing model: %s \n', label_drug{idrug})

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

num_unit_plot = length(unique(att_table.unit(idx_unit)));

fprintf('\t Found %d units \n', num_unit_plot)

for idp = 1:length(data2plot)

% plot violinplot

iplot = iplot+1;

% init axes

subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

if idrug==1

axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);

else

axH = plotj_initAx(fSet);

end

hold on

% extract smaller table

% v_att_table = unstack(att_table(idx_unit,:), data2plot{idp}, {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {data2plot{idp},'drug_on','attention','unit','attention_drug_on'}));

switch data2plot{idp}

case {'gain_log', 'gain'}

gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit'}, 'VariableNamingRule', 'preserve');

variable_names = gv.Properties.VariableNames;

for iatt_drug = 1:length(labels_violin)

mus = [gv.(labels_violin{iatt_drug}).mus];

s2s = [gv.(labels_violin{iatt_drug}).s2s];

mus = mus(:);

s2s = s2s(:);

% get binomial fit

rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);

gv_pop = GainVariability;

gv_pop.rhat = rhat;

m2v = gv_pop.pred_variance(xval_gain);

% plot

scatter(mus, s2s, fSet.MarkerSize, colors_violin(iatt_drug,:), 'filled', 'MarkerFaceAlpha', 0.5);

plot(m2v.mean, m2v.variance,'Color', colors_violin(iatt_drug,:), 'linew', fSet.LineWidth);

xlim([0.8 max(xval_gain)])

ylim([0.8 max(xval_gain)])

plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)

set(gca,'xscale','log','yscale','log')

end

% axis square

% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)

if idrug==length(label_drug)

xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)

end

ylabel('Variance (spikes^{2})','FontSize',fSet.Fontsize_text)

axis square

iplot = iplot+1;

subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

if idrug==1

axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);

else

axH = plotj_initAx(fSet);

end

hold on

end

% plot violin

v = violinplot(v_att_table.(data2plot{idp}), v_att_table.attention_drug_on, 'ShowMean', true);

XTickLabel = get(gca,'XTickLabel');

for iv = 1:length(v)

v(iv).ScatterPlot.SizeData = 10;

idx = strcmpi(XTickLabel(iv), labels_violin);

v(iv).ViolinColor = colors_violin(idx,:);

end

% plot lines for each unit

[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));

for iv = 1:length(v)

% data_x(:,iv) = v(iv).ScatterPlot.XData;

% data_y(:,iv) = v(iv).ScatterPlot.YData;

if ~any(isnan(v_att_table.(data2plot{idp})))

% works if data is complete

data_x(:,iv) = v(iv).ScatterPlot.XData;

data_y(:,iv) = v(iv).ScatterPlot.YData;

else

% catch if nan values are present

idx = strcmpi(XTickLabel(iv), labels_violin);

idx_trial = v_att_table.attention_drug_on==labels_violin{idx};

data_x(:,iv) = iv;

data_y(:,iv) = v_att_table.(data2plot{idp})(idx_trial);

end

plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);

% fit model

lme_population = fitlme(v_att_table,[data2plot{idp} ' ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding

disp(lme_population)

coefficients = lme_population.Coefficients;

betas = coefficients.Estimate(2:end);

betas_se = coefficients.SE(2:end);

p_val = coefficients.pValue(2:end);

% p_string = get_significance_strings(p_val, 'rounded', 0);

p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});

p_string{end+1} = sprintf('n=%d',num_unit_plot);

% % full betastring: Estimate, SE, p-value

% betastring = {...

% sprintf('drug: %s = %1.2f%s%1.2f, %s', '\beta', betas(1), '\pm', betas_se(1), p_string{1}), ...

% sprintf('att: %s = %1.2f%s%1.2f, %s', '\beta', betas(2), '\pm', betas_se(2), p_string{2}), ...

% sprintf('att %s drug: %s = %1.2f%s%1.2f, %s', '\times', '\beta', betas(3), '\pm', betas_se(3), p_string{3}), ...

% sprintf('n = %d',length(unique(v_att_table.unit)))};

% % part betastring: Estimate, SE, p-value

% betastring = {...

% sprintf('drug: %s = %1.2f, %s', '\beta', betas(1), p_string{1}), ...

% sprintf('att: %s = %1.2f, %s', '\beta', betas(2), p_string{2}), ...

% sprintf('att %s drug: %s = %1.2f, %s', '\times', '\beta', betas(3), p_string{3}), ...

% sprintf('n = %d',length(unique(v_att_table.unit)))};

% plot text

YLIM = get(gca,'ylim');

switch data2plot{idp}

case {'gain_log', 'gain'}

set(gca, 'ylim', YLIM .* [1 1.5])

otherwise

set(gca, 'ylim', YLIM .* [1 1.15])

end

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.95);

y_pos = get_value_range(tmp_y, 0.95);

h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right', 'FontSize', fSet.Fontsize_text);

[p_masked] = p_val<0.05;

p_masked = [p_masked ; false];

plotj_text_emphasise(h_text, 1, 'italic', p_masked);

plotj_text_emphasise(h_text, 1, 'bold', p_masked);

switch data2plot{idp}

case 'rate'

if idp==1

ylabel({label_drug{idrug}, 'Firing rate (Hz)'}, 'FontSize', fSet.Fontsize_text)

else

ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)

end

case 'FF'

ylabel('Fano Factor', 'FontSize', fSet.Fontsize_text)

case 'gain'

ylabel('\sigma_{G}^{2}', 'FontSize', fSet.Fontsize_text)

case 'gain_log'

ylabel('log(\sigma_{G}^{2})', 'FontSize', fSet.Fontsize_text)

end

if idrug == 2

set(gca,'XTickLabelRotation',25)

else

set(gca, 'XTickLabel', '')

end

% attention/drug auroc

iplot = iplot+1;

% init axes

subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

if idrug==1

axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);

else

axH = plotj_initAx(fSet);

end

hold on

% get data - AUROC

idx_cond = att_table.attention=='Attend RF';

idx = idx_unit & idx_cond;

tmp_data = zeros(length(find(idx))/2,1);

tmp_data(:,1) = att_table.roc(idx & att_table.drug_on=='Drug off');

tmp_data(:,2) = att_table.roc(idx & att_table.drug_on=='Drug on');

plotj_scatter(tmp_data, ...

'markerStyle', markerstyle, 'MarkerSize', markersize, ...

'markerFaceColor', [0.5 0.5 0.5], 'markerFaceAlpha', 0.5, ...

'markerEdgeColor', [0.5 0.5 0.5], 'markerEdgeAlpha', 0.5, ...

'axislimit', [0 1]);

% stats

P(idrug) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');

axis square

p_string = get_significance_strings(P(idrug), 'rounded', 0);

% mean difference

delta_data = tmp_data(:,2)-tmp_data(:,1);

% effect size

d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');

% print result

fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

'AUROC', ...

label_drug{idrug}, ...

'AUROC', ...

mean(delta_data), ...

std(delta_data)/sqrt(size(tmp_data,1)), ...

p_string, ...

% set text

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.05);

y_pos = get_value_range(tmp_y, 0.95);

h_text(idrug) = text(x_pos, y_pos, ...

sprintf('%s (n=%d)', p_string, size(tmp_data,1)), ...

'FontSize', fSet.Fontsize_text);

xlabel('Attention AUROC no drug','FontSize', fSet.Fontsize_text)

ylabel('Attention AUROC drug','FontSize', fSet.Fontsize_text)

end

Testing model: Dopamine 

	 Found 31 units 

Linear mixed-effects model fit by ML

Model information:
    Number of observations             124
    Fixed effects coefficients           4
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    rate ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    848.12    865.04    -418.06          836.12  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate    SE         tStat       DF     pValue        Lower      Upper  
    {'(Intercept)'                        }         19.132      2.7502      6.9564    120    1.9744e-10     13.687     24.577
    {'drug_on_Drug on'                    }        -2.3639     0.38775     -6.0964    120    1.3585e-08    -3.1316    -1.5962
    {'attention_Attend RF'                }         2.6939     0.38775      6.9476    120    2.0643e-10     1.9262     3.4616
    {'drug_on_Drug on:attention_Attend RF'}        -0.2884     0.38775    -0.74379    120       0.45846    -1.0561    0.47931

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower    Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        15.16       11.76    19.543

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        4.3178      3.7398    4.9851

Linear mixed-effects model fit by ML

Model information:
    Number of observations             124
    Fixed effects coefficients           4
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    FF ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    447.27    464.19    -217.63          435.27  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate    SE         tStat       DF     pValue        Lower       Upper  
    {'(Intercept)'                        }          3.4901    0.21273      16.406    120    1.9163e-32      3.0689     3.9112
    {'drug_on_Drug on'                    }        -0.02305    0.10547    -0.21855    120       0.82737    -0.23186    0.18576
    {'attention_Attend RF'                }         0.10797    0.10547      1.0238    120       0.30801    -0.10084    0.31678
    {'drug_on_Drug on:attention_Attend RF'}        0.093561    0.10547     0.88713    120       0.37678    -0.11525    0.30238

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.0286      0.73701    1.4355

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        1.1744      1.0172    1.3559

Linear mixed-effects model fit by ML

Model information:
    Number of observations             118
    Fixed effects coefficients           4
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    gain_log ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    238.62    255.24    -113.31          226.62  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate     SE          tStat      DF     pValue        Lower        Upper     
    {'(Intercept)'                        }          -1.8048     0.14524    -12.426    114    6.1881e-23      -2.0925       -1.5171
    {'drug_on_Drug on'                    }          0.19251    0.042312     4.5498    114    1.3492e-05      0.10869       0.27633
    {'attention_Attend RF'                }        -0.088546    0.042312    -2.0927    114      0.038593     -0.17236    -0.0047272
    {'drug_on_Drug on:attention_Attend RF'}         0.023579    0.042629    0.55313    114       0.58126    -0.060868       0.10803

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.77318     0.58935    1.0144

Group: Error
    Name               Estimate    Lower      Upper  
    {'Res Std'}        0.4571      0.39407    0.53022

AUROC: Dopamine, delta-AUROC -0.02 +- 0.01, p = 0.094, Cohens's d = -0.28

Testing model: SCH23390 

	 Found 14 units 

Linear mixed-effects model fit by ML

Model information:
    Number of observations              56
    Fixed effects coefficients           4
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    rate ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    355.56    367.71    -171.78          343.56  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate    SE         tStat       DF    pValue        Lower      Upper   
    {'(Intercept)'                        }          15.922     2.6476       6.014    52    1.8299e-07      10.61      21.235
    {'drug_on_Drug on'                    }         -1.3651    0.44484     -3.0688    52     0.0034106    -2.2578    -0.47248
    {'attention_Attend RF'                }          2.5624    0.44484      5.7602    52    4.5931e-07     1.6697       3.455
    {'drug_on_Drug on:attention_Attend RF'}        -0.16397    0.44484    -0.36861    52       0.71391    -1.0566     0.72866

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        9.7654      6.6701    14.297

Group: Error
    Name               Estimate    Lower    Upper 
    {'Res Std'}        3.3289      2.688    4.1226

Linear mixed-effects model fit by ML

Model information:
    Number of observations              56
    Fixed effects coefficients           4
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    FF ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    232.04    244.19    -110.02          220.04  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate    SE         tStat     DF    pValue        Lower        Upper  
    {'(Intercept)'                        }         3.8125     0.39513    9.6487    52    3.5217e-13       3.0196     4.6054
    {'drug_on_Drug on'                    }         0.2818     0.19273    1.4622    52       0.14971     -0.10494    0.66855
    {'attention_Attend RF'                }        0.30423     0.19273    1.5785    52       0.12051    -0.082515    0.69097
    {'drug_on_Drug on:attention_Attend RF'}        0.21862     0.19273    1.1343    52       0.26185     -0.16812    0.60536

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.2907      0.79022    2.108

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        1.4423      1.1646    1.7862

Linear mixed-effects model fit by ML

Model information:
    Number of observations              53
    Fixed effects coefficients           4
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    gain_log ~ 1 + drug_on*attention + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    131.34    143.16    -59.668          119.34  

Fixed effects coefficients (95% CIs):
    Name                                           Estimate    SE          tStat      DF    pValue        Lower        Upper   
    {'(Intercept)'                        }         -1.5478     0.28729    -5.3876    49    2.0228e-06      -2.1251    -0.97044
    {'drug_on_Drug on'                    }         0.24421    0.071712     3.4054    49     0.0013266       0.1001     0.38832
    {'attention_Attend RF'                }        -0.14014    0.071712    -1.9542    49      0.056402     -0.28425    0.003974
    {'drug_on_Drug on:attention_Attend RF'}         0.15056     0.07154     2.1046    49      0.040483    0.0067948     0.29433

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.0411      0.70072    1.5467

Group: Error
    Name               Estimate    Lower      Upper  
    {'Res Std'}        0.51583     0.41309    0.64411

AUROC: SCH23390, delta-AUROC -0.02 +- 0.03, p = 0.583, Cohens's d = -0.20

Post hoc comparisons: scatter plots and signed-rank tests

plotting and post-hoc tests

ncol = 2 * (length(label_drug));

nrow = 3; %

[fH, fSet] = plotj_initFig('figNum', 2, 'width', 20, 'height', 12, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap .* [1 0.8];

[P, h_text] = deal( NaN(length(label_drug), nrow, 2, 2) ); % drugtype, actvity-type, attention, drug-offon,

iplot=0;

for irow = 1:3

if irow == 1

% mean rate

datatype = 'rate';

ylim2use = [0.7 100];

scale ='Log';

x_text = 0.005;

y_text = [0.9];

elseif irow == 2

% Fano factor

datatype = 'FF';

ylim2use = [0 10];

scale ='Linear';

x_text = 0.1;

y_text = [0.95];

elseif irow == 3

% gain variability

datatype = 'gain_log';

ylim2use = [-4 2];

scale ='Linear';

x_text = 0.1;

y_text = [0.9];

end

data2plot = att_table.(datatype);

for idrug = 1:length(label_drug)

for iatt = 1:2

iplot = iplot+1;

subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

axH = plotj_initAx(fSet, 'axlabelDisplacement', [0.05, 0.02]);

hold on

if irow==1 && iatt==1

% sup-title

set(gca,'units','normalized')

text(100, 400, ['\bf ' label_drug{idrug}], 'FontSize', fSet.Fontsize_title, 'Interpreter', 'tex')

end

if (irow==1)

title(sprintf('%s', label_attention{iatt}), 'FontSize', fSet.Fontsize_title)

end

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_att = att_table.attention == label_attention{iatt};

unit2plot = idx_unit ...

& idx_att;

% extract data

tmp_data = [data2plot(unit2plot & att_table.drug_on==label_drug_onoff{1}) ...

data2plot(unit2plot & att_table.drug_on==label_drug_onoff{2})];

plotj_scatter(tmp_data, ...

'MarkerSize', markersize, ...

'markerFaceColor', fSet.colors(1,:), 'markerFaceAlpha', 0.5, ...

'axislimit', ylim2use);

% stats

P(idrug,irow,iatt) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');

% effect size

d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');

% mean difference

delta_data = tmp_data(:,2)-tmp_data(:,1);

axis square

p_string = get_significance_strings(P(idrug,irow,iatt), 'rounded', 0);

% print result

fprintf('%s: %s, attend %s: delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

datatype, ...

label_drug{idrug}, ...

label_attention{iatt}, ...

datatype, ...

nanmean(delta_data), ...

nanstd(delta_data)/sqrt(length(tmp_data)), ...

p_string, ...

% plot text

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, x_text);

y_pos = get_value_range(tmp_y, y_text);

h_text(idrug,irow,iatt) = text(x_pos, y_pos, ...

sprintf('%s', p_string), ...

'Color', fSet.colors(1,:));

if iatt==2

x_pos = get_value_range(tmp_x, 0.6);

y_pos = get_value_range(tmp_y, 0.2);

text(x_pos, y_pos, sprintf('(n=%d)', length(tmp_data)), 'Color', fSet.colors(1,:));

end

xlim(ylim2use)

ylim(ylim2use)

set(gca,'YScale', scale, 'XScale', scale)

% % set(gca,'', [1 10 100])

if irow == 1

set(gca,'Xtick', [1 10 100], 'XTickLabel', [1 10 100])

set(gca,'Ytick', [1 10 100], 'YTickLabel', [1 10 100])

xlabel('Firing rate no drug (Hz)','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Firing rate drug (Hz)','FontSize', fSet.Fontsize_text)

end

elseif irow == 2

xlabel('Fano Factor no drug','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Fano Factor drug','FontSize', fSet.Fontsize_text)

end

elseif irow == 3

xlabel('Gain variability no drug','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Gain variability drug','FontSize', fSet.Fontsize_text)

end

rate: Dopamine, attend Attend RF: delta-rate 5.30 +- 0.87, p < 0.001, Cohens's d = 1.09
rate: Dopamine, attend Attend away: delta-rate 4.15 +- 0.84, p < 0.001, Cohens's d = 0.88
rate: SCH23390, attend Attend RF: delta-rate 3.06 +- 0.80, p = 0.002, Cohens's d = 1.03
rate: SCH23390, attend Attend away: delta-rate 2.40 +- 0.91, p = 0.017, Cohens's d = 0.71
FF: Dopamine, attend Attend RF: delta-FF -0.14 +- 0.36, p = 0.779, Cohens's d = -0.07
FF: Dopamine, attend Attend away: delta-FF 0.23 +- 0.21, p = 0.504, Cohens's d = 0.20
FF: SCH23390, attend Attend RF: delta-FF -1.00 +- 0.79, p = 0.296, Cohens's d = -0.34
FF: SCH23390, attend Attend away: delta-FF -0.13 +- 0.33, p = 0.903, Cohens's d = -0.10
gain_log: Dopamine, attend Attend RF: delta-gain_log -0.46 +- 0.14, p = 0.006, Cohens's d = -0.57
gain_log: Dopamine, attend Attend away: delta-gain_log -0.34 +- 0.11, p = 0.006, Cohens's d = -0.59
gain_log: SCH23390, attend Attend RF: delta-gain_log -0.82 +- 0.31, p = 0.009, Cohens's d = -0.82
gain_log: SCH23390, attend Attend away: delta-gain_log -0.18 +- 0.12, p = 0.244, Cohens's d = -0.15

[p_fdr, p_masked] = FDR(P, 0.05);

plotj_text_emphasise(h_text, p_masked, 'italic');

plotj_text_emphasise(h_text, p_masked, 'bold');

plot histogram cell-widths

colors = get_colors('spikewidth');

binsize = 25;

binedges = 50:binsize:600;

assert(any(ismember(binedges,waveform_cutoff)), 'threshold needs to be included in binedges')

[fH, fSet] = plotj_initFig('width', 20, 'height', 12, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap .* [1 0.8];

ncol=4;

nrow=2;

subtightplot(nrow, ncol, 1, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

plotj_initAx(fSet, 'axlabel', 1, 'axlabelDisplacement', [0.07, 0.01]);

hold on

x = waveform_time-(7*6*5.4);

plot(x, unitlist.waveform(unitlist.unit_class=='Narrow',:)*-1, 'Color', [colors(1,:) 0.5])

plot(x, unitlist.waveform(unitlist.unit_class=='Broad',:)*-1, 'Color', [colors(2,:) 0.5])

xlabel(['Time from peak of spike (' plotj_symbol('mu') 's)'], 'FontSize', fSet.Fontsize_text)

ylabel(['Normalized voltage'], 'FontSize', fSet.Fontsize_text)

xlim(minmax(waveform_time)+[0 -200]-7*6*5.4)

axis square

subtightplot(nrow, ncol, 2, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

plotj_initAx(fSet, 'axlabel', 2, 'axlabelDisplacement', [0.07, 0.01]);

hold on

h = plotj_hist({unitlist.peak_to_trough_time(unitlist.unit_class=='Narrow'), unitlist.peak_to_trough_time(unitlist.unit_class=='Broad')}, ...

'bins', binedges, ...

'EdgeColor', colors, ...

'histstyle', 'stairs',...

'plotmean', 0, ...

'LineWidth', fSet.LineWidth);

xlabel(['Peak to through time (' char(181) 's)'], 'FontSize', fSet.Fontsize_text)

ylabel('Number of cells', 'FontSize', fSet.Fontsize_text)

axis square

hleg = legend(h, {'Narrow','Broad'}, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'AutoUpdate', 'Off', 'Location', 'East');

hleg.Position = hleg.Position + [0.08 0 0 0];

% do hartigan's dip test to see whether there is a significant dip in this

% distribution

% [dip, p_value, xlow,xup]=HartigansDipSignifTest(allSpikeWaveform.width(idx_unit),5000);

[dip, p_value, xlow, xup, boot_dip]=CalibratedHartigansDipSignifTest(unitlist.peak_to_trough_time, 10000);

p_string = get_significance_strings(p_value, 'round', 0);

h_text(1) = text(400, 13, {'Calibrated Hartigan''s dip test:', p_string}, 'FontSize', fSet.Fontsize_text);

LME model fit

model fit to 'rate', 'FF' and 'gain' across attention and drug conditions, for broad and narrow spiking cells. We use random intercepts for each unit to account for the repeated measures

% Note these excellent explanations of effect coded predictors:
% - https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-interpret-the-coefficients-of-an-effect-coded-variable-involved-in-an-interaction-in-a-regression-model/#effect-coded-predictors-interacted-with-a-continuous-covariate
% - https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faqwhat-is-effect-coding/
data2plot = {'rate', 'FF', 'gain_log'};
for idrug = 1:num_drug
    
    idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
    
    tmp_table = att_table(idx_unit,:);
    
    for idp = 1:length(data2plot)
        
        fprintf('\nDRUG: %s, DATATYPE: %s, UNIT SELECTION: %s, n=%d \n', label_drug{idrug}, data2plot{idp}, selectivity_criterium, length(unique(att_table.unit(idx_unit))))
        
        % model fit across attention, drug and unit_class
        mod_population.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention*drug_on*unit_class + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
        mod_population.anova = anova(mod_population.lme);
        
        %         % hierarchical model fits
        %         mod.m1_intercept.lme = fitlme(tmp_table,[data2plot{idp} ' ~ 1 + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m2_att.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m3_drug.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m4_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m5_att_drug.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention*drug_on + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m6_att_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m7_drug_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + drug_on:unit_class + (1|unit)'],'DummyVarCoding','effects');
        %         mod.m8_att_drug_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + drug_on:unit_class + attention:drug_on:unit_class + (1|unit)'],'DummyVarCoding','effects');
        %
        %         % likelihood ratio tests of hierarchical model fits
        %         mod.m2_att.compare = compare(mod.m1_intercept.lme, mod.m2_att.lme, 'CheckNesting', true);
        %         mod.m3_drug.compare = compare(mod.m2_att.lme, mod.m3_drug.lme, 'CheckNesting', true);
        %         mod.m4_unitc.compare = compare(mod.m3_drug.lme, mod.m4_unitc.lme, 'CheckNesting', true);
        %         mod.m5_att_drug.compare = compare(mod.m4_unitc.lme, mod.m5_att_drug.lme, 'CheckNesting', true);
        %         mod.m6_att_unitc.compare = compare(mod.m5_att_drug.lme, mod.m6_att_unitc.lme, 'CheckNesting', true);
        %         mod.m7_drug_unitc.compare = compare(mod.m6_att_unitc.lme, mod.m7_drug_unitc.lme, 'CheckNesting', true);
        %         mod.m8_att_drug_unitc.compare = compare(mod.m7_drug_unitc.lme, mod.m8_att_drug_unitc.lme, 'CheckNesting', true);
        %
        % print stats for each group
        disp(grpstats(tmp_table,{'attention','drug_on','unit_class'},{'min','max','mean','median'},'DataVars',data2plot{idp}))
        
        % print coefficients
        disp(mod_population.lme)
    end
end
DRUG: Dopamine, DATATYPE: rate, UNIT SELECTION: att&dru, n=31 
                                    attention     drug_on     unit_class    GroupCount    min_rate    max_rate    mean_rate    median_rate
                                   ___________    ________    __________    __________    ________    ________    _________    ___________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow          14        0.79125      47.567      16.942        15.312   
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           17         3.4929      69.675      21.011        10.703   
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow          14         1.5655      52.707      20.546        21.285   
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           17         4.5532      80.411      27.717        15.746   
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow          14         1.9689      34.203      12.246        9.0068   
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           17         4.1932       54.01      16.105        8.9194   
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow          14         4.2369      38.811      15.463        13.631   
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           17          5.637       60.66      21.026        11.448   
Linear mixed-effects model fit by ML

Model information:
    Number of observations             124
    Fixed effects coefficients           8
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    rate ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    852.27    880.47    -416.13          832.27  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate    SE         tStat       DF     pValue        Lower       Upper  
    {'(Intercept)'                                          }          18.882      2.724      6.9317    116    2.4955e-10      13.487     24.277
    {'drug_on_Drug on'                                      }         -2.3057    0.38341     -6.0137    116    2.1581e-08     -3.0651    -1.5463
    {'attention_Attend RF'                                  }          2.6719    0.38341      6.9688    116    2.0728e-10      1.9125     3.4313
    {'unit_class_Narrow'                                    }         -2.5828      2.724    -0.94818    116       0.34501     -7.9781     2.8124
    {'drug_on_Drug on:attention_Attend RF'                  }         -0.2715    0.38341     -0.7081    116        0.4803     -1.0309     0.4879
    {'drug_on_Drug on:unit_class_Narrow'                    }         0.60072    0.38341      1.5668    116       0.11989    -0.15868     1.3601
    {'attention_Attend RF:unit_class_Narrow'                }        -0.22728    0.38341    -0.59279    116       0.55448    -0.98668    0.53212
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}          0.1747    0.38341     0.45564    116        0.6495     -0.5847     0.9341

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        14.945      11.593    19.266

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        4.2495      3.6806    4.9062
DRUG: Dopamine, DATATYPE: FF, UNIT SELECTION: att&dru, n=31 
                                    attention     drug_on     unit_class    GroupCount    min_FF    max_FF    mean_FF    median_FF
                                   ___________    ________    __________    __________    ______    ______    _______    _________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow          14        2.0911    7.1985    3.8655      3.0195  
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           17        1.2696    6.3605    3.5064      3.3962  
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow          14        1.4285    5.2381    2.8399      2.3452  
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           17        2.0925    10.152    4.0938      3.2168  
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow          14        1.7752    5.4315     3.213      2.8758  
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           17        2.0443    5.7579    3.3087      2.6984  
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow          14        1.8161     6.673    3.5405      3.2446  
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           17        1.4416     8.299    3.4643      2.8067  
Linear mixed-effects model fit by ML

Model information:
    Number of observations             124
    Fixed effects coefficients           8
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    FF ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    446.19    474.39    -213.09          426.19  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate      SE         tStat        DF     pValue       Lower        Upper   
    {'(Intercept)'                                          }             3.479    0.21274       16.353    116    6.766e-32       3.0576      3.9004
    {'drug_on_Drug on'                                      }        -0.0056168    0.10107    -0.055573    116      0.95578      -0.2058     0.19457
    {'attention_Attend RF'                                  }          0.097384    0.10107      0.96352    116      0.33729      -0.1028     0.29757
    {'unit_class_Narrow'                                    }          -0.11428    0.21274     -0.53717    116      0.59218     -0.53564     0.30708
    {'drug_on_Drug on:attention_Attend RF'                  }           0.11515    0.10107       1.1393    116      0.25692    -0.085032     0.31533
    {'drug_on_Drug on:unit_class_Narrow'                    }           0.18014    0.10107       1.7823    116     0.077315    -0.020044     0.38032
    {'attention_Attend RF:unit_class_Narrow'                }           -0.1094    0.10107      -1.0824    116      0.28132     -0.30958    0.090784
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}           0.22309    0.10107       2.2073    116     0.029263     0.022908     0.42327

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.0374      0.75014    1.4346

Group: Error
    Name               Estimate    Lower      Upper 
    {'Res Std'}        1.1202      0.97024    1.2933
DRUG: Dopamine, DATATYPE: gain_log, UNIT SELECTION: att&dru, n=31 
                                    attention     drug_on     unit_class    GroupCount    min_gain_log    max_gain_log    mean_gain_log    median_gain_log
                                   ___________    ________    __________    __________    ____________    ____________    _____________    _______________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow          14          -3.1554            1.0697        -1.4737           -1.4938    
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           17          -4.2117          -0.46886        -1.8043           -1.9479    
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow          14          -3.2009        -0.0074051         -2.115           -2.3693    
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           17          -3.6859          -0.82216        -2.0801            -2.031    
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow          14          -2.9678           0.86773        -1.2674           -1.4614    
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           17          -3.4723          -0.53745        -1.6883           -1.6183    
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow          14          -3.2715          -0.11799        -1.6461           -1.8262    
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           17          -3.5353          -0.76188        -2.0821           -2.0195    
Linear mixed-effects model fit by ML

Model information:
    Number of observations             118
    Fixed effects coefficients           8
    Random effects coefficients         31
    Covariance parameters                2

Formula:
    gain_log ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    240.21    267.92    -110.1           220.21  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate     SE          tStat      DF     pValue        Lower        Upper    
    {'(Intercept)'                                          }          -1.7965     0.14453     -12.43    110    1.0898e-22      -2.0829      -1.5101
    {'drug_on_Drug on'                                      }          0.20057      0.0414     4.8446    110    4.1833e-06      0.11852      0.28261
    {'attention_Attend RF'                                  }        -0.099717      0.0414    -2.4086    110      0.017675     -0.18176    -0.017672
    {'unit_class_Narrow'                                    }          0.11883     0.14453    0.82221    110       0.41274     -0.16759      0.40526
    {'drug_on_Drug on:attention_Attend RF'                  }         0.033293    0.041782    0.79683    110       0.42726    -0.049509      0.11609
    {'drug_on_Drug on:unit_class_Narrow'                    }         0.034814      0.0414     0.8409    110       0.40223    -0.047232      0.11686
    {'attention_Attend RF:unit_class_Narrow'                }         -0.06958      0.0414    -1.6807    110      0.095665     -0.15163     0.012465
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}          0.06441    0.041782     1.5416    110       0.12605    -0.018392      0.14721

Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.76695     0.58502    1.0055

Group: Error
    Name               Estimate    Lower      Upper  
    {'Res Std'}        0.44218     0.38117    0.51295
DRUG: SCH23390, DATATYPE: rate, UNIT SELECTION: att&dru, n=14 
                                    attention     drug_on     unit_class    GroupCount    min_rate    max_rate    mean_rate    median_rate
                                   ___________    ________    __________    __________    ________    ________    _________    ___________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow           3         10.527      38.003      19.858        11.043   
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           11         1.8333      31.818      16.164        16.361   
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow           3         12.542      43.356      23.049        13.248   
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           11         2.7802        37.4      19.186        17.697   
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow           3         4.9403      21.982      11.365        7.1728   
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           11         1.0897       26.58      12.375        8.2434   
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow           3         6.5644      24.866      12.802         6.976   
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           11         1.6288      34.826      15.041        11.932   
Linear mixed-effects model fit by ML

Model information:
    Number of observations              56
    Fixed effects coefficients           8
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    rate ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    356.62    376.87    -168.31          336.62  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate    SE         tStat       DF    pValue        Lower       Upper   
    {'(Intercept)'                                          }           16.23     3.2229      5.0358    48    7.1363e-06      9.7498       22.71
    {'drug_on_Drug on'                                      }         -1.2895    0.49925     -2.5828    48      0.012903     -2.2933    -0.28567
    {'attention_Attend RF'                                  }          3.3342    0.49925      6.6784    48    2.2901e-08      2.3304       4.338
    {'unit_class_Narrow'                                    }          0.5384     3.2229     0.16705    48       0.86803     -5.9418      7.0186
    {'drug_on_Drug on:attention_Attend RF'                  }        -0.26378    0.49925    -0.52836    48       0.59968     -1.2676     0.74002
    {'drug_on_Drug on:unit_class_Narrow'                    }         0.13238    0.49925     0.26515    48       0.79203    -0.87143      1.1362
    {'attention_Attend RF:unit_class_Narrow'                }          1.3506    0.49925      2.7054    48     0.0094162     0.34684      2.3545
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}        -0.17466    0.49925    -0.34986    48       0.72798     -1.1785     0.82914

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        9.7769      6.6891    14.29

Group: Error
    Name               Estimate    Lower     Upper
    {'Res Std'}        3.066       2.4757    3.797
DRUG: SCH23390, DATATYPE: FF, UNIT SELECTION: att&dru, n=14 
                                    attention     drug_on     unit_class    GroupCount    min_FF    max_FF    mean_FF    median_FF
                                   ___________    ________    __________    __________    ______    ______    _______    _________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow           3        4.8438    6.9749    5.6723      5.1983  
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           11        1.2078    14.527    4.3294      3.6886  
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow           3        4.7386    7.6174    5.9872      5.6056  
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           11        1.2523    4.7912    2.9697      2.8743  
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow           3        3.2601     4.786    4.1195      4.3125  
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           11        2.1601    7.0603     3.422      2.7141  
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow           3        4.0334    5.6038    4.7933      4.7427  
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           11        1.8187    4.2039    3.0774      3.1568  
Linear mixed-effects model fit by ML

Model information:
    Number of observations              56
    Fixed effects coefficients           8
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    FF ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    233.16    253.42    -106.58          213.16  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate     SE         tStat       DF    pValue        Lower         Upper  
    {'(Intercept)'                                          }           4.2964    0.42499      10.109    48    1.7726e-13        3.4419     5.1509
    {'drug_on_Drug on'                                      }         0.089451    0.22558     0.39654    48       0.69347      -0.36411    0.54301
    {'attention_Attend RF'                                  }           0.4433    0.22558      1.9652    48      0.055197     -0.010258    0.89686
    {'unit_class_Narrow'                                    }          0.84672    0.42499      1.9923    48       0.05204    -0.0077845     1.7012
    {'drug_on_Drug on:attention_Attend RF'                  }          0.17175    0.22558     0.76135    48       0.45017      -0.28181    0.62531
    {'drug_on_Drug on:unit_class_Narrow'                    }         -0.33662    0.22558     -1.4922    48       0.14218      -0.79018    0.11694
    {'attention_Attend RF:unit_class_Narrow'                }          0.24338    0.22558      1.0789    48       0.28602      -0.21018    0.69694
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}        -0.082033    0.22558    -0.36365    48       0.71771      -0.53559    0.37153

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.106       0.65592    1.8649

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        1.3853      1.1186    1.7156
DRUG: SCH23390, DATATYPE: gain_log, UNIT SELECTION: att&dru, n=14 
                                    attention     drug_on     unit_class    GroupCount    min_gain_log    max_gain_log    mean_gain_log    median_gain_log
                                   ___________    ________    __________    __________    ____________    ____________    _____________    _______________

    Attend RF_Drug on_Narrow       Attend RF      Drug on       Narrow           3          -2.0016         -0.65222         -1.2007          -0.94844    
    Attend RF_Drug on_Broad        Attend RF      Drug on       Broad           11          -3.5623           1.1763          -1.121           -1.6503    
    Attend RF_Drug off_Narrow      Attend RF      Drug off      Narrow           3          -1.9344         -0.89988         -1.3305           -1.1571    
    Attend RF_Drug off_Broad       Attend RF      Drug off      Broad           11          -4.0364         -0.22889         -2.2878           -2.6556    
    Attend away_Drug on_Narrow     Attend away    Drug on       Narrow           3          -1.9583         -0.54467         -1.0418          -0.62238    
    Attend away_Drug on_Broad      Attend away    Drug on       Broad           11          -3.1992           1.2952         -1.3882            -1.477    
    Attend away_Drug off_Narrow    Attend away    Drug off      Narrow           3           -1.679         -0.38222        -0.92504          -0.71391    
    Attend away_Drug off_Broad     Attend away    Drug off      Broad           11          -3.2254          0.82684         -1.5237           -1.3357    
Linear mixed-effects model fit by ML

Model information:
    Number of observations              53
    Fixed effects coefficients           8
    Random effects coefficients         14
    Covariance parameters                2

Formula:
    gain_log ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    134.64    154.34    -57.318          114.64  

Fixed effects coefficients (95% CIs):
    Name                                                             Estimate      SE          tStat        DF    pValue        Lower         Upper    
    {'(Intercept)'                                          }           -1.3918      0.3414      -4.0767    45    0.00018358       -2.0794     -0.70419
    {'drug_on_Drug on'                                      }           0.16056    0.080691       1.9898    45      0.052716    -0.0019651      0.32308
    {'attention_Attend RF'                                  }          -0.13652    0.080691      -1.6918    45       0.09759      -0.29904     0.026004
    {'unit_class_Narrow'                                    }            0.2673      0.3414      0.78294    45       0.43776      -0.42033      0.95492
    {'drug_on_Drug on:attention_Attend RF'                  }           0.12026    0.080615       1.4918    45       0.14272     -0.042104      0.28263
    {'drug_on_Drug on:unit_class_Narrow'                    }          -0.15731    0.080691      -1.9496    45      0.057475      -0.31983    0.0052072
    {'attention_Attend RF:unit_class_Narrow'                }        -0.0045813    0.080691    -0.056775    45       0.95498       -0.1671      0.15794
    {'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'}         -0.058643    0.080615     -0.72744    45       0.47073      -0.22101      0.10372

Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.0187      0.68687    1.5108

Group: Error
    Name               Estimate    Lower      Upper  
    {'Res Std'}        0.48979     0.39227    0.61156

Violin, population plot: rate, FF and gain, resolved for unit class

data2plot = {'rate','FF', 'gain_log'};

[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');

colors_violin = reshape(colors_violin, [4,3]);

labels_violin = reshape(labels_violin, [4,1]);

labels_violin_short = reshape(labels_violin_short, [4,1]);

% plotting and post-hoc tests

ncol = num_drug*length(label_unitclass); % data2plot x unit

nrow = length(data2plot);

iplot = 0;

% use subplot indices to combine violin and post-hoc scatter plots

idx_axlabel = [1 NaN nrow+1 NaN 2 NaN nrow+2 NaN 3 NaN nrow+3 NaN];

[fH, fSet] = plotj_initFig('figNum', 3, 'width', 20, 'height', 15, 'Journal', plot_conventions);

fSet.subplotGap = [0.1, 0.06];

for idp = 1:length(data2plot)

for idrug = 1:num_drug

% clear param

axpos = zeros(length(label_unitclass),4);

% fit model

fprintf('Plot data for drug: %s \n', label_drug{idrug})

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

fprintf('\t Found %d units \n', length(unique(att_table.unit(idx_unit))))

% plot, violinplot for each unitclass

for iunitc = 1:length(label_unitclass) % for each unit_class

iplot = iplot+1;

% init axes

subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

if iunitc==1

axH = plotj_initAx(fSet, 'axlabel', idx_axlabel(iplot), 'axlabelDisplacement', [0.05, 0.02]);

else

axH = plotj_initAx(fSet);

axH.Position = axH.Position - [0.08 0 0 0];

end

hold on

h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);

% find units in class

idx_unit_class = att_table.unit_class==label_unitclass(iunitc);

num_unit_plot = length(unique(att_table.unit(idx_unit & idx_unit_class)));

% extract smaller table

v_att_table = att_table(idx_unit & idx_unit_class, ismember(att_table.Properties.VariableNames, {data2plot{idp},'drug_on','attention','unit','attention_drug_on'}));

fprintf('\t Found %d independent units in unit_class: %s \n', length(unique(v_att_table.unit)), label_unitclass{iunitc})

% plot violin

v = violinplot(v_att_table.(data2plot{idp}), v_att_table.attention_drug_on, 'ShowMean', true);

XTickLabel = get(gca,'XTickLabel');

for iv = 1:length(v)

v(iv).ScatterPlot.SizeData = 10;

idx = strcmpi(XTickLabel(iv), labels_violin);

v(iv).ViolinColor = colors_violin(idx,:);

end

% plot lines for each unit

[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));

for iv = 1:length(v)

if ~any(isnan(v_att_table.(data2plot{idp})))

% works if data is complete

data_x(:,iv) = v(iv).ScatterPlot.XData;

data_y(:,iv) = v(iv).ScatterPlot.YData;

else

% catch if nan values are present

idx = strcmpi(XTickLabel(iv), labels_violin);

idx_trial = v_att_table.attention_drug_on==labels_violin{idx};

data_x(:,iv) = iv;

data_y(:,iv) = v_att_table.(data2plot{idp})(idx_trial);

end

plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);

% fit model

lme_rate = fitlme(v_att_table,[data2plot{idp} ' ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding

% fprintf('Model fit: %s, unit: %s\n', label_drug{idrug}, label_unittype{iunitc})

% disp(lme_rate)

coefficients = lme_rate.Coefficients;

betas = coefficients.Estimate(2:end);

betas_se = coefficients.SE(2:end);

p_val = coefficients.pValue(2:end);

% p_string = get_significance_strings(p_val, 'rounded', 0);

p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});

p_string{end+1} = sprintf('n=%d',num_unit_plot);

% % full betastring: Estimate, SE, p-value

% betastring = {...

% sprintf('drug: %s = %1.2f%s%1.2f, %s', '\beta', betas(1), '\pm', betas_se(1), p_string{1}), ...

% sprintf('att: %s = %1.2f%s%1.2f, %s', '\beta', betas(2), '\pm', betas_se(2), p_string{2}), ...

% sprintf('att %s drug: %s = %1.2f%s%1.2f, %s', '\times', '\beta', betas(3), '\pm', betas_se(3), p_string{3}), ...

% sprintf('n = %d',length(unique(v_att_table.unit)))};

% % part betastring: Estimate, SE, p-value

% betastring = {...

% sprintf('drug: %s = %1.2f, %s', '\beta', betas(1), p_string{1}), ...

% sprintf('att: %s = %1.2f, %s', '\beta', betas(2), p_string{2}), ...

% sprintf('att %s drug: %s = %1.2f, %s', '\times', '\beta', betas(3), p_string{3}), ...

% sprintf('n = %d',length(unique(v_att_table.unit)))};

% plot text

YLIM = get(gca,'ylim');

switch data2plot{idp}

case 'rate'

set(gca, 'ylim', [0 1000])

y_text = 0.2;

case 'FF'

set(gca, 'ylim', YLIM .* [1 1.5])

y_text = 0.8;

case 'gain_log'

set(gca, 'ylim', [-5 5])

y_text = 0.8;

end

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 1);

y_pos = get_value_range(tmp_y, y_text);

h_text = text(x_pos, y_pos, p_string, ...

'HorizontalAlignment','right',...

'FontSize',fSet.Fontsize_text_in);

if iunitc==1

switch data2plot{idp}

case 'rate'

ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)

case 'FF'

ylabel('Fano Factor', 'FontSize', fSet.Fontsize_text)

case 'gain'

ylabel('\sigma_{G}^{2}', 'FontSize', fSet.Fontsize_text)

case 'gain_log'

ylabel('log(\sigma_{G}^{2})', 'FontSize', fSet.Fontsize_text)

end

if (idp<length(data2plot)) || (idrug>1) || (iunitc>1)

set(gca,'XTickLabel',[])

else

set(gca,'XTickLabelRotation',-25)

end

% [p_fdr, p_masked] = FDR(p_val, 0.05);

[p_masked] = p_val<0.05;

p_masked = [p_masked ; false];

plotj_text_emphasise(h_text, 1, 'italic', p_masked);

plotj_text_emphasise(h_text, 1, 'bold', p_masked);

switch data2plot{idp}

case {'rate','gain'}

set(gca,'YScale','log')

end

Plot data for drug: Dopamine 

	 Found 31 units 

	 Found 14 independent units in unit_class: Narrow 
	 Found 17 independent units in unit_class: Broad 

Plot data for drug: SCH23390 

	 Found 14 units 

	 Found 3 independent units in unit_class: Narrow 
	 Found 11 independent units in unit_class: Broad 

Plot data for drug: Dopamine 

	 Found 31 units 

	 Found 14 independent units in unit_class: Narrow 
	 Found 17 independent units in unit_class: Broad 

Plot data for drug: SCH23390 

	 Found 14 units 

	 Found 3 independent units in unit_class: Narrow 
	 Found 11 independent units in unit_class: Broad 

Plot data for drug: Dopamine 

	 Found 31 units 

	 Found 14 independent units in unit_class: Narrow 
	 Found 17 independent units in unit_class: Broad 

Plot data for drug: SCH23390 

	 Found 14 units 

	 Found 3 independent units in unit_class: Narrow 
	 Found 11 independent units in unit_class: Broad 

% create suptitle

ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized','clipping','off');

xpos_idx = 1:2:ncol;

for idrug = 1:num_drug

xpos = (1/(ncol)) * xpos_idx(idrug);

ht = text(xpos, 0.95, [label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title * 1.2, 'HorizontalAlignment', 'center');

ht.FontWeight = 'bold';

end

Post hoc comparisons: scatter plots and signed-rank tests, resolved for unit type

plotting and post-hoc tests

ncol = 2 * (length(label_drug));

nrow = 3; %

markerstyle = {'o','v'};

colors_scatter = get_colors('spikewidth');

idx_subplot = [1 2 ncol+1 ncol+2 2*ncol+1 2*ncol+2; 3 4 ncol+3 ncol+4 2*ncol+3 2*ncol+4] ;

idx_axlabel = [1 NaN 2 NaN 3 NaN; 4 NaN 5 NaN 6 NaN];

[fH, fSet] = plotj_initFig('figNum', 2, 'width', 20, 'height', 16, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap .* [1 0.8];

[P, h_text] = deal( NaN(length(label_drug), nrow, 2, 2) ); % drugtype, actvity-type, attention, drug-offon,

for idrug = 1:length(label_drug)

iplot=0;

for irow = 1:nrow

if irow == 1

% mean rate

datatype = 'rate';

ylim2use = [0.7 120];

scale ='Log';

x_text = 0.005;

y_text = [0.9 0.4];

elseif irow == 2

% Fano factor

datatype = 'FF';

ylim2use = [0 4];

scale ='Linear';

x_text = 0.1;

y_text = [0.95 0.1];

elseif irow == 3

% Gain variability

datatype = 'gain_log';

ylim2use = [-4 2];

scale ='Linear';

x_text = 0.1;

y_text = [0.95 0.1];

end

data2plot = att_table.(datatype);

for iatt = 1:2

iplot = iplot+1;

subtightplot(nrow, ncol, idx_subplot(idrug, iplot), fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

axH = plotj_initAx(fSet, 'axlabel', idx_axlabel(idrug, iplot), 'axlabelDisplacement', [0.05, 0.02]);

hold on

if irow==1 && iatt==1

% sup-title

set(gca,'units','normalized')

text(100, 400, ['\bf ' label_drug{idrug}], 'FontSize', fSet.Fontsize_title, 'Interpreter', 'tex')

end

title(sprintf('Attend %s', label_attention{iatt}), 'FontSize', fSet.Fontsize_title)

for iunitc = 1:2

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_att = att_table.attention == label_attention{iatt};

idx_unit_class = att_table.unit_class == label_unitclass{iunitc};

unit2plot = idx_unit ...

& idx_att ...

& idx_unit_class;

% extract data

tmp_data = [data2plot(unit2plot & att_table.drug_on==label_drug_onoff{1}) ...

data2plot(unit2plot & att_table.drug_on==label_drug_onoff{2})];

plotj_scatter(tmp_data, ...

'markerStyle', markerstyle(iunitc), 'MarkerSize', markersize, ...

'markerFaceColor', colors_scatter(iunitc,:), 'markerFaceAlpha', 0.5, ...

'markerEdgeColor', colors_scatter(iunitc,:), 'markerEdgeAlpha', 0.5, ...

'axislimit', ylim2use);

% stats

P(idrug,irow,iatt,iunitc) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');

% effect size

d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');

% mean difference

delta_data = tmp_data(:,2)-tmp_data(:,1);

axis square

p_string = get_significance_strings(P(idrug,irow,iatt,iunitc), 'rounded', 0);

% print result

fprintf('%s: %s, %s, attend %s: delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

datatype, ...

label_drug{idrug}, ...

label_unitclass{iunitc}, ...

label_attention{iatt}, ...

datatype, ...

mean(delta_data), ...

std(delta_data)/sqrt(length(tmp_data)), ...

p_string, ...

% plot text

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, x_text);

y_pos = get_value_range(tmp_y, y_text(1)-y_text(2)*(iunitc-1));

h_text(idrug,irow,iatt,iunitc) = text(x_pos, y_pos, ...

sprintf('%s', p_string), ...

'Color', colors_scatter(iunitc,:));

if iatt==2

x_pos = get_value_range(tmp_x, 0.6);

y_pos = get_value_range(tmp_y, 0.1+0.1*(iunitc-1));

text(x_pos, y_pos, sprintf('%s (n=%d)', label_unitclass{iunitc}, length(tmp_data)), 'Color', colors_scatter(iunitc,:));

end

xlim(ylim2use)

ylim(ylim2use)

set(gca,'YScale', scale, 'XScale', scale)

% % set(gca,'', [1 10 100])

if strcmpi(datatype, 'rate')

set(gca,'Xtick', [1 10 100], 'XTickLabel', [1 10 100])

set(gca,'Ytick', [1 10 100], 'YTickLabel', [1 10 100])

xlabel('Firing rate no drug (Hz)','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Firing rate drug (Hz)','FontSize', fSet.Fontsize_text)

end

elseif strcmp(datatype, 'FF')

xlabel('Fano Factor no drug','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Fano Factor drug','FontSize', fSet.Fontsize_text)

end

elseif strcmp(datatype, 'gain_log')

xlabel('Gain no drug','FontSize', fSet.Fontsize_text)

if idrug==1 && iatt==1

ylabel('Gain drug','FontSize', fSet.Fontsize_text)

end

rate: Dopamine, Narrow, attend Attend RF: delta-rate 3.60 +- 1.12, p = 0.009, Cohens's d = 0.86
rate: Dopamine, Broad, attend Attend RF: delta-rate 6.71 +- 1.22, p < 0.001, Cohens's d = 1.33
rate: Dopamine, Narrow, attend Attend away: delta-rate 3.22 +- 1.06, p = 0.009, Cohens's d = 0.81
rate: Dopamine, Broad, attend Attend away: delta-rate 4.92 +- 1.26, p < 0.001, Cohens's d = 0.95
FF: Dopamine, Narrow, attend Attend RF: delta-FF -1.03 +- 0.58, p = 0.119, Cohens's d = -0.47
FF: Dopamine, Broad, attend Attend RF: delta-FF 0.59 +- 0.39, p = 0.174, Cohens's d = 0.37
FF: Dopamine, Narrow, attend Attend away: delta-FF 0.33 +- 0.30, p = 0.326, Cohens's d = 0.30
FF: Dopamine, Broad, attend Attend away: delta-FF 0.16 +- 0.31, p = 0.927, Cohens's d = 0.12
gain_log: Dopamine, Narrow, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.032, Cohens's d = -0.64
gain_log: Dopamine, Broad, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.093, Cohens's d = -0.54
gain_log: Dopamine, Narrow, attend Attend away: delta-gain_log NaN +- NaN, p = 0.204, Cohens's d = -0.55
gain_log: Dopamine, Broad, attend Attend away: delta-gain_log -0.39 +- 0.15, p = 0.011, Cohens's d = -0.65
rate: SCH23390, Narrow, attend Attend RF: delta-rate 3.19 +- 1.14, p = 0.250, Cohens's d = 1.62
rate: SCH23390, Broad, attend Attend RF: delta-rate 3.02 +- 0.99, p = 0.010, Cohens's d = 0.92
rate: SCH23390, Narrow, attend Attend away: delta-rate 1.44 +- 0.89, p = 0.500, Cohens's d = 0.93
rate: SCH23390, Broad, attend Attend away: delta-rate 2.67 +- 1.13, p = 0.042, Cohens's d = 0.71
FF: SCH23390, Narrow, attend Attend RF: delta-FF 0.31 +- 0.39, p = 0.500, Cohens's d = 0.47
FF: SCH23390, Broad, attend Attend RF: delta-FF -1.36 +- 0.98, p = 0.206, Cohens's d = -0.42
FF: SCH23390, Narrow, attend Attend away: delta-FF 0.67 +- 0.39, p = 0.500, Cohens's d = 1.00
FF: SCH23390, Broad, attend Attend away: delta-FF -0.34 +- 0.39, p = 0.638, Cohens's d = -0.26
gain_log: SCH23390, Narrow, attend Attend RF: delta-gain_log -0.13 +- 0.10, p = 0.500, Cohens's d = -0.76
gain_log: SCH23390, Broad, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.020, Cohens's d = -0.94
gain_log: SCH23390, Narrow, attend Attend away: delta-gain_log 0.12 +- 0.11, p = 0.500, Cohens's d = 0.62
gain_log: SCH23390, Broad, attend Attend away: delta-gain_log NaN +- NaN, p = 0.160, Cohens's d = -0.27

[p_fdr, p_masked] = FDR(P, 0.05);

plotj_text_emphasise(h_text, p_masked, 'italic');

plotj_text_emphasise(h_text, p_masked, 'bold');

plot gain, mean to variance relationship, separate for unit-types

[colors, labels] = get_colors('att_drug');

[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');

colors_violin = reshape(colors_violin, [4,3]);

labels_violin = reshape(labels_violin, [4,1]);

ncol = length(label_drug);

nrow = 2; % unit types

[fH, fSet] = plotj_initFig('width',20,'height',15,'Journal',plot_conventions);

clear legendText h

gain_table_full = [];

iplot=0;

xval_gain = 1:500;

for idrug = 1:length(label_drug)

for iunitc = 1:length(unique(unitlist.unit_class))

iplot = iplot+1;

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);

% subplot axis

h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

if iunitc==1

plotj_initAx(fSet, 'axlabel', idrug);

else

plotj_initAx(fSet);

end

hold on

if idrug==1

h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);

end

gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

variable_names = gv.Properties.VariableNames;

for iatt_drug = 1:length(labels_violin)

mus = [gv.(labels_violin{iatt_drug}).mus];

s2s = [gv.(labels_violin{iatt_drug}).s2s];

mus = mus(:);

s2s = s2s(:);

% get binomial fit

rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);

gv_pop = GainVariability;

gv_pop.rhat = rhat;

m2v = gv_pop.pred_variance(xval_gain);

% plot

h(iatt) = scatter(mus,s2s,fSet.MarkerSize*3,colors_violin(iatt_drug,:),'filled','MarkerFaceAlpha',0.5);

plot(m2v.mean,m2v.variance,'Color',colors_violin(iatt_drug,:),'linew',fSet.LineWidth);

xlim([0.8 max(xval_gain)])

ylim([0.8 max(xval_gain)])

plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)

set(gca,'xscale','log','yscale','log')

end

% axis square

% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)

if idrug==length(label_drug)

xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)

end

if iunitc==1

ylabel({[label_drug{idrug} label_drug_ext{idrug}],'Variance (spikes^{2})'},'FontSize',fSet.Fontsize_text)

end

% violin

% init axes

pos = get(h_ax, 'Position');

h_ax_inset = axes('Position', [pos(1)+0.22 pos(2)+0.03 0.16 0.15]) ; % inset

plotj_initAx(fSet);

hold on

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);

num_unit_plot = length(unique(att_table.unit(idx_unit)));

fprintf('\t Found %d units \n', num_unit_plot)

% extract smaller table

v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {'gain_log','drug_on','attention','unit','attention_drug_on'}));

% plot violin

v = violinplot(v_att_table.gain_log, v_att_table.attention_drug_on, 'ShowMean', true);

XTickLabel = get(gca,'XTickLabel');

for iv = 1:length(v)

v(iv).ScatterPlot.SizeData = 10;

idx = strcmpi(XTickLabel(iv), labels_violin);

v(iv).ViolinColor = colors_violin(idx,:);

end

set(gca,'XTickLabel',[])

h_ax_inset.XAxis.Visible = 'off';

% plot lines for each unit

[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));

for iv = 1:length(v)

if ~any(isnan(v_att_table.gain_log))

% works if data is complete

data_x(:,iv) = v(iv).ScatterPlot.XData;

data_y(:,iv) = v(iv).ScatterPlot.YData;

else

% catch if nan values are present

idx = strcmpi(XTickLabel(iv), labels_violin);

idx_trial = v_att_table.attention_drug_on==labels_violin{idx};

data_x(:,iv) = iv;

data_y(:,iv) = v_att_table.gain_log(idx_trial);

end

plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);

% fit model

lme_population = fitlme(v_att_table,['gain_log ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding

disp(lme_population.Coefficients)

coefficients = lme_population.Coefficients;

betas = coefficients.Estimate(2:end);

betas_se = coefficients.SE(2:end);

p_val = coefficients.pValue(2:end);

% p_string = get_significance_strings(p_val, 'rounded', 0);

p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});

p_string{end+1} = ['n = ' num2str(num_unit_plot)];

% plot text

YLIM = get(gca,'ylim');

set(gca, 'ylim', YLIM .* [1 1.25])

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.95);

y_pos = get_value_range(tmp_y, 1.3);

h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right');

ylabel('log(\sigma^{2}_G )', 'FontSize', fSet.Fontsize_text)

% [p_fdr, p_masked] = FDR(p_val, 0.05);

[p_masked] = p_val<0.05;

p_masked = [p_masked ; false];

plotj_text_emphasise(h_text, 1, 'italic', p_masked);

plotj_text_emphasise(h_text, 1, 'bold', p_masked);

end

	 Found 14 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate    SE          tStat      DF    pValue        Lower        Upper    
    {'(Intercept)'                        }         -1.6769     0.24109    -6.9553    47    9.5535e-09      -2.1619      -1.1919
    {'drug_on_Drug on'                    }         0.23563    0.074521     3.1619    47      0.002744     0.085713      0.38555
    {'attention_Attend RF'                }        -0.16931    0.074521     -2.272    47      0.027706     -0.31923    -0.019394
    {'drug_on_Drug on:attention_Attend RF'}        0.097133    0.075694     1.2832    47       0.20571    -0.055145      0.24941

	 Found 17 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate     SE          tStat       DF    pValue        Lower       Upper   
    {'(Intercept)'                        }          -1.9154      0.1718     -11.149    63    1.5163e-16     -2.2587      -1.572
    {'drug_on_Drug on'                    }          0.16573    0.045026      3.6809    63    0.00048392    0.075757     0.25571
    {'attention_Attend RF'                }        -0.030159    0.045026    -0.66981    63       0.50543    -0.12014    0.059818
    {'drug_on_Drug on:attention_Attend RF'}        -0.031138    0.045026    -0.69156    63       0.49176    -0.12111    0.058839

	 Found 3 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate     SE          tStat      DF    pValue       Lower         Upper    
    {'(Intercept)'                        }          -1.1245     0.31895    -3.5257    8     0.0077812         -1.86     -0.38901
    {'drug_on_Drug on'                    }        0.0032421    0.029576    0.10962    8       0.91541     -0.064961     0.071446
    {'attention_Attend RF'                }          -0.1411    0.029576    -4.7706    8     0.0014075       -0.2093    -0.072895
    {'drug_on_Drug on:attention_Attend RF'}          0.06162    0.029576     2.0834    8      0.070737    -0.0065831      0.12982

	 Found 11 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate    SE          tStat      DF    pValue        Lower        Upper   
    {'(Intercept)'                        }        -1.6587      0.34668    -4.7845    37    2.7399e-05      -2.3611    -0.95625
    {'drug_on_Drug on'                    }        0.31791     0.088241     3.6027    37    0.00092093      0.13911      0.4967
    {'attention_Attend RF'                }        -0.1319     0.088241    -1.4948    37       0.14346     -0.31069    0.046895
    {'drug_on_Drug on:attention_Attend RF'}        0.17933     0.087884     2.0405    37      0.048482    0.0012596      0.3574

Attend ROC, drug/no drug, narrow/broad. DrugMI-ejecCurrent

% datatype for dose-response curve

datatype = 'MI';

% datatype = 'AUROC';

% datatype = 'gain_log';

markerstyle = {'o','o'};

ncol = length(label_drug) ;

nrow = 1;

[fH, fSet] = plotj_initFig('width', 15, 'height', 8, 'Journal',plot_conventions);

switch datatype

case 'MI'

v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);

iunitc=[];

ylabel2use = 'Drug modulation index';

ylim2use = {[-0.5 0.5],[-0.4 0.2]};

case 'AUROC'

% get data - AUROC

idx_cond = att_table.attention=='Attend RF';

tmp_data = zeros(length(find(idx_cond))/2,1);

tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');

tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');

data2plot = diff(tmp_data,1,2);

iunitc = [];

ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];

ylim2use = {[-0.25 0.3],[-0.25 0.3]};

case 'gain_log'

v_att_table = unstack(att_table, 'gain_log', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

% difference in attentional modulation of gain variabiliy

data2plot1 = (v_att_table.("Attend RF, Drug on") - v_att_table.("Attend away, Drug on")) ./ (v_att_table.("Attend RF, Drug on") + v_att_table.("Attend away, Drug on"));

data2plot2 = (v_att_table.("Attend RF, Drug off") - v_att_table.("Attend away, Drug off")) ./ (v_att_table.("Attend RF, Drug off") + v_att_table.("Attend away, Drug off"));

data2plot = data2plot1 - data2plot2;

iunitc=[];

end

clear ax h

for idrug = 1:length(label_drug)

clear text_legend

idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});

ax(idrug) = subtightplot(nrow, ncol, idrug, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

plotj_initAx(fSet, 'axlabel', idrug);

hold on

title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

% ejectionCurrent

xlabel_event = 'Ejection Current';

xlim2use = [17 90];

tmp_data = [unitlist.EjectCurrent(idx_unit) data2plot(idx_unit)];

hscat = plotj_scatter(tmp_data, ...

'markerStyle', markerstyle, 'markerSize', markersize, ...

'markerFaceColor', 'k', 'markerFaceAlpha', 0.5, ...

'axislimit', [], ...

'unityLine', 0);

text_legend{1} = ['(n=' num2str(length(find(idx_unit))) ')'];

data_table = table(unitlist.EjectCurrent_centered(idx_unit), data2plot(idx_unit), 'VariableNames', {'x', 'y'});

[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);

predict_data.x = predict_data.x + mean(unique(unitlist.EjectCurrent));

switch datatype

case 'MI'

filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criterium));

case 'AUROC'

filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criterium));

case 'gain_log'

end

writetable(data_table, filename)

% fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)

if any(strcmp('fit_mean', predict_data.Properties.VariableNames))

h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);

cond1 = exist([filename(1:end-4) '_R.csv'], 'file');

cond2 = strcmpi(model,'quadratic');

if cond1 && cond2

utable = readtable([filename(1:end-4) '_R.csv']);

if utable.usig==1

pstring_chi = [pstring_chi ' (U^{+})'];

else

pstring_chi = [pstring_chi ' (U^{-})'];

end

else

pstring_chi = '';

end

switch model

case 'linear'

betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];

% model_color = colors(2,:);

case 'quadratic'

betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];

otherwise

betaString = [];

end

axis square

xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)

if idrug==1

ylabel(ylabel2use, 'FontSize', fSet.Fontsize_text)

end

if ~isempty(ylim2use)

ylim(ylim2use{idrug})

end

if ~isempty(xlim2use)

xlim(xlim2use)

end

% set ax param

set(gca,'YDir','reverse')

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.05);

y_pos = get_value_range(tmp_y, 0.15);

text(x_pos, y_pos, {betaString, pstring_chi, text_legend{1}}, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3])

end

Attend ROC, drug/no drug, resolved for narrow/broad. DrugMI-ejecCurrent

% datatype for dose-response curve

datatype = 'MI';

% datatype = 'AUROC';

% datatype = 'gain_log';

colors = get_colors('spikewidth');

markerstyle = {'o','o'};

idx_group = 1;

ncol = length(label_drug);

nrow = 2;

[fH, fSet] = plotj_initFig('width', 20, 'height', 14, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap.*[1.2 .8];

axislim = [0.25 0.9];

[P, h_text] = deal( NaN(length(label_drug), 2) ); % drugtype, unit-type,

[P_roc, h_text_roc] = deal( NaN(length(label_drug), 2) ); % drugtype, unit-type,

iplot=0;

for idrug = 1:length(label_drug)

iplot = iplot+1;

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

% subplot axis

h_ax_subplot = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

plotj_initAx(fSet, 'axlabel', idrug, 'axlabelDisplacement', [0.01, 0.02]);

hold on

title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

% inset axis

pos = get(h_ax_subplot, 'Position');

h_ax_inset = axes('Position', [pos(1)+0.28 pos(2)+0.035 0.1 0.11]) ; % inset

plotj_initAx(fSet);

hold on

% lme

for iunitc = 1:2

set(fH, 'currentaxes', h_ax_subplot);

unit2plot = idx_unit & att_table.unit_class==label_unitclass(iunitc);

% get data - AUROC

idx_cond = att_table.attention=='Attend RF';

idx = unit2plot & idx_cond;

tmp_data = zeros(length(find(idx))/2,1);

tmp_data(:,1) = att_table.roc(idx & att_table.drug_on=='Drug off');

tmp_data(:,2) = att_table.roc(idx & att_table.drug_on=='Drug on');

% [bf10,p,CI,stats] = bf.ttest(tmp_data(:,1),tmp_data(:,2))

plotj_scatter(tmp_data, ...

'markerStyle', markerstyle, 'MarkerSize', markersize, ...

'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...

'markerEdgeColor', colors(iunitc,:), 'markerEdgeAlpha', 0.5, ...

'axislimit', axislim);

% stats

P(idrug,iunitc) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');

axis square

p_string = get_significance_strings(P(idrug,iunitc), 'rounded', 0);

% mean difference

delta_data = tmp_data(:,2)-tmp_data(:,1);

% effect size

d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');

% print result

fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

'AUROC', ...

label_drug{idrug}, ...

'AUROC', ...

mean(delta_data), ...

std(delta_data)/sqrt(size(tmp_data,1)), ...

p_string, ...

% set text

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.05);

y_pos = get_value_range(tmp_y, 0.95 - 0.07 * (iunitc-1));

h_text(idrug,iunitc) = text(x_pos, y_pos, ...

sprintf('%s (n=%d)', p_string, size(tmp_data,1)), ...

'Color', colors(iunitc,:), 'FontSize', fSet.Fontsize_ax);

% inset

% -----

idx = tmp_data(:,1)<0.5;

tmp_data(idx,:) = 1-tmp_data(idx,:);

% [bf10,p,CI,stats] = bf.ttest(tmp_data(:,1),tmp_data(:,2))

P_roc(idrug,iunitc) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');

bar_data = diff(tmp_data,1,2);

bar_y = mean(bar_data);

bar_y_err = std(bar_data)/sqrt(length(bar_data));

% plot

set(fH, 'currentaxes', h_ax_inset);

h_scatter = scatter(iunitc + (rand(length(bar_data),1)-0.5)*0.1, bar_data, ...

'.', 'MarkerFaceColor', colors(iunitc,:),'MarkerEdgeColor', colors(iunitc,:));

h_scatter.MarkerFaceAlpha = 0.5;

h_scatter.MarkerEdgeAlpha = 0.5;

plot(iunitc + [-0.25 0.25], [bar_y bar_y], 'Color', colors(iunitc,:), 'linew', fSet.LineWidth+1)

plot(iunitc + [0 0], bar_y + [-bar_y_err bar_y_err], 'Color', colors(iunitc,:), 'linew', fSet.LineWidth)

p_string = get_significance_strings(P_roc(idrug,iunitc), 'rounded', 0);

% mean difference

delta_data = tmp_data(:,2)-tmp_data(:,1);

% effect size

d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');

% print result

fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

'AUROC-CORRECT', ...

label_drug{idrug}, ...

'AUROC', ...

mean(delta_data), ...

std(delta_data)/sqrt(size(tmp_data,1)), ...

p_string, ...

% set text

% tmp_y = get(gca,'ylim');

%y_pos = get_value_range(tmp_y, 0.95);

y_pos = 0.1*idrug;

h_text_roc(idrug,iunitc) = text(iunitc+0.1, y_pos, ...

p_string, ...

'Color', colors(iunitc,:), ...

'FontSize', fSet.Fontsize_ax);

end

set(h_ax_inset, 'XTick', [1 2])

% ylabel(h_ax_inset, [plotj_symbol('Delta') ' | AUROC-0.5 |'],'FontSize', fSet.Fontsize_text)

ylabel(h_ax_inset, [plotj_symbol('Delta') ' AUROC-c'],'FontSize', fSet.Fontsize_text)

xlabel(h_ax_subplot, 'Attention AUROC no drug','FontSize', fSet.Fontsize_text)

if idrug==1

ylabel(h_ax_subplot, 'Attention AUROC drug','FontSize', fSet.Fontsize_text)

end

AUROC: Dopamine, delta-AUROC 0.00 +- 0.02, p = 0.952, Cohens's d = 0.03

AUROC-CORRECT: Dopamine, delta-AUROC -0.01 +- 0.02, p = 0.542, Cohens's d = -0.16

AUROC: Dopamine, delta-AUROC -0.03 +- 0.01, p = 0.009, Cohens's d = -0.69

AUROC-CORRECT: Dopamine, delta-AUROC -0.02 +- 0.01, p = 0.080, Cohens's d = -0.38

AUROC: SCH23390, delta-AUROC -0.03 +- 0.04, p = 0.750, Cohens's d = -0.54

AUROC-CORRECT: SCH23390, delta-AUROC -0.03 +- 0.04, p = 0.750, Cohens's d = -0.54

AUROC: SCH23390, delta-AUROC -0.02 +- 0.04, p = 0.700, Cohens's d = -0.15

AUROC-CORRECT: SCH23390, delta-AUROC -0.04 +- 0.03, p = 0.278, Cohens's d = -0.33

% FDR correction

[p_fdr, p_masked] = FDR(P, 0.05);

[p_fdr, p_masked_corr] = FDR(P_roc, 0.05);

plotj_text_emphasise(h_text, p_masked, 'italic');

plotj_text_emphasise(h_text, p_masked, 'bold');

% drugMI/attAUROC against ejection current

switch datatype

case 'MI'

v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);

iunitc=[];

ylabel2use = 'Drug modulation index';

ylim2use = {[-0.5 0.5],[-0.4 0.2]};

case 'AUROC'

% get data - AUROC

idx_cond = att_table.attention=='Attend RF';

tmp_data = zeros(length(find(idx_cond))/2,1);

tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');

tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');

data2plot = diff(tmp_data,1,2);

iunitc = [];

ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];

ylim2use = {[-0.25 0.3],[-0.25 0.3]};

case 'gain_log'

v_att_table = unstack(att_table, 'gain_log', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

% difference in attentional modulation of gain variabiliy

data2plot1 = (v_att_table.("Attend RF, Drug on") - v_att_table.("Attend away, Drug on")) ./ (v_att_table.("Attend RF, Drug on") + v_att_table.("Attend away, Drug on"));

data2plot2 = (v_att_table.("Attend RF, Drug off") - v_att_table.("Attend away, Drug off")) ./ (v_att_table.("Attend RF, Drug off") + v_att_table.("Attend away, Drug off"));

data2plot = data2plot1 - data2plot2;

iunitc=[];

end

clear ax h

for idrug = 1:ncol

clear text_legend

iplot=iplot+1;

idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});

ax(idrug) = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

plotj_initAx(fSet, 'axlabel', iplot, 'axlabelDisplacement', [0.01, 0.02]);

hold on

title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

if ~isempty(iunitc)

unit2plot = idx_unit & unitlist.unit_class==label_unitclass(iunitc);

else

unit2plot = idx_unit;

end

length(find(unit2plot))

% ejectionCurrent

xlabel_event = 'Ejection Current';

xlim2use = [17 90];

if length(unique(unitlist.unit_class(unit2plot)))==1

tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];

hscat = plotj_scatter(tmp_data, ...

'markerStyle', markerstyle, 'markerSize', markersize, ...

'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...

'axislimit', [], ...

'unityLine', 0);

text_legend{1} = [label_unitclass{iunitc} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==iunitc))) ')'];

else

tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];

hscat = plotj_scatter(tmp_data, ...

'dataIndex', grp2idx(unitlist.unit_class(unit2plot)), ...

'markerStyle', markerstyle, 'MarkerSize', markersize, ...

'markerFaceColor', (colors(1:2,:)), 'markerFaceAlpha', [0.5 0.5], ...

'markerEdgeColor', (colors(1:2,:)), 'markerEdgeAlpha', [0.5 0.5], ...

'axislimit', [], ...

'unityLine', 0);

text_legend{1} = [label_unitclass{1} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(1)))) ')'];

text_legend{2} = [label_unitclass{2} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(2)))) ')'];

end

data_table = table(unitlist.EjectCurrent_centered(unit2plot), data2plot(unit2plot), 'VariableNames', {'x', 'y'});

[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);

predict_data.x = predict_data.x + mean(unique(unitlist.EjectCurrent));

switch datatype

case 'MI'

filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criterium));

case 'AUROC'

filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criterium));

case 'gain_log'

end

writetable(data_table, filename)

% fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)

if any(strcmp('fit_mean', predict_data.Properties.VariableNames))

h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);

cond1 = exist([filename(1:end-4) '_R.csv'], 'file');

cond2 = strcmpi(model,'quadratic');

if cond1 && cond2

utable = readtable([filename(1:end-4) '_R.csv']);

if utable.usig==1

pstring_chi = [pstring_chi ' (U^{+})'];

else

pstring_chi = [pstring_chi ' (U^{-})'];

end

else

pstring_chi = '';

end

switch model

case 'linear'

betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];

% model_color = colors(2,:);

case 'quadratic'

betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];

otherwise

betaString = [];

end

axis square

xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)

if idrug==1

ylabel(ylabel2use, 'FontSize', fSet.Fontsize_text)

end

if ~isempty(ylim2use)

ylim(ylim2use{idrug})

end

if ~isempty(xlim2use)

xlim(xlim2use)

end

% plot text indicating number of units

set(gca, 'Units', 'normalized');

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.6);

y_pos(1) = get_value_range(tmp_y, 0.8);

y_pos(2) = get_value_range(tmp_y, 0.9);

text(x_pos, y_pos(1), text_legend{1}, 'FontSize', fSet.Fontsize_text, 'Color', colors(1,:))

text(x_pos, y_pos(2), text_legend{2}, 'FontSize', fSet.Fontsize_text, 'Color', colors(2,:))

% set ax param

set(gca,'YDir','reverse')

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.05);

y_pos = get_value_range(tmp_y, 0.1);

text(x_pos, y_pos, {betaString, pstring_chi }, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3])

end

ans = 31

ans = 14

DrugMI-ejecCurrent across selectivity indices

colors = get_colors('spikewidth');

datatype = 'MI';

datatype = 'AUROC';

selectivity_criteria = {'none', 'att', 'dru', 'att&dru'};

label_criteria = {'No subselection', 'Attention-selective', 'Drug-selective', {'Attention & Drug', 'selective'}};

ncol = length(selectivity_criteria);

nrow = length(label_drug);

iplot = 0;

[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap.*[1 .75];

idx_subplot = [1:ncol ; (ncol+1):ncol*2];

% drugMI/attAUROC against ejection current

switch datatype

case 'MI'

v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);

iunitc=[];

ylabel2use = 'Drug modulation index';

ylim2use = [];

case 'AUROC'

% data2plot = squeeze(diff(rate_ROC.roc,1,2));

% data2plot = squeeze(mean(data2plot(:,idx_attention_roc==1),2)); % average over relevant auroc conditions

idx_cond = att_table.attention=='Attend RF';

tmp_data = zeros(length(find(idx_cond))/2,1);

tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');

tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');

data2plot = diff(tmp_data,1,2);

iunitc = [];

ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];

ylim2use = [-0.25 0.3];

end

[P,h_text] = deal(zeros(length(selectivity_criteria), length(label_drug)));

for icrit = 1:length(selectivity_criteria)

selectivity_criterium = selectivity_criteria{icrit};

clear ax h

for idrug = 1:length(label_drug)

clear text_legend

iplot=iplot+1;

idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});

ax(idrug) = subtightplot(nrow, ncol, idx_subplot(iplot), fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

if idrug==1

plotj_initAx(fSet, 'axlabel', icrit);

else

plotj_initAx(fSet);

end

hold on

if idrug==1

% title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

h_title = title(label_criteria{icrit}, 'FontSize', fSet.Fontsize_title);

set(h_title, 'Units', 'normalized');

h_title.Position = h_title.Position+[0 0.2 0];

end

if ~isempty(iunitc)

unit2plot = idx_unit & unitlist.unit_class==label_unitclass{iunitc};

else

unit2plot = idx_unit;

end

% ejectionCurrent

xlabel_event = 'Ejection Current';

xlim2use = [17 95];

if length(unique(unitlist.unit_class(unit2plot)))==1

tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];

hscat = plotj_scatter(tmp_data, ...

'markerStyle', {'.'}, 'markerSize', markersize, ...

'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...

'axislimit', [], ...

'unityLine', 0);

text_legend{1} = [label_unitclass{iunitc} ' (n=' num2str(length(find(unit_class(idx_unit)==iunitc))) ')'];

else

tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];

hscat = plotj_scatter(tmp_data, ...

'dataIndex', grp2idx(unitlist.unit_class(unit2plot)), ...

'markerStyle', {'o','o'}, 'MarkerSize', markersize, ...

'markerFaceColor', (colors(1:2,:)), 'markerFaceAlpha', [0.5 0.5], ...

'markerEdgeColor', (colors(1:2,:)), 'markerEdgeAlpha', [0.5 0.5], ...

'axislimit', [], ...

'unityLine', 0);

text_legend{1} = [label_unitclass{1} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(1)))) ')'];

text_legend{2} = [label_unitclass{2} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(2)))) ')'];

end

data_table = table(unitlist.EjectCurrent(unit2plot), data2plot(unit2plot), 'VariableNames', {'x', 'y'});

[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);

P(icrit,idrug) = stats2report.pValue;

switch datatype

case 'MI'

filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criteria{icrit}));

case 'AUROC'

filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criteria{icrit}));

end

writetable(data_table, filename);

fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)

if any(strcmp('fit_mean', predict_data.Properties.VariableNames))

h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);

[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);

cond1 = exist([filename(1:end-4) '_R.csv'], 'file');

cond2 = strcmpi(model,'quadratic');

if cond1 && cond2

utable = readtable([filename(1:end-4) '_R.csv']);

if utable.usig==1

pstring_chi = [pstring_chi ' (U^{+})'];

else

pstring_chi = [pstring_chi ' (U^{-})'];

end

else

pstring_chi = '';

end

switch model

case 'linear'

betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];

% model_color = colors(2,:);

case 'quadratic'

betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];

otherwise

betaString = [];

end

axis square

if idrug==length(label_drug)

xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)

end

if icrit==1

ylabel({label_drug{idrug},ylabel2use}, 'FontSize', fSet.Fontsize_text)

end

if ~isempty(ylim2use)

ylim(ylim2use)

end

if ~isempty(xlim2use)

xlim(xlim2use)

end

% plot text indicating number of units

set(gca, 'Units', 'normalized');

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.6);

y_pos(1) = get_value_range(tmp_y, 0.8);

y_pos(2) = get_value_range(tmp_y, 0.9);

text(x_pos, y_pos(1), text_legend{1}, 'FontSize', fSet.Fontsize_text, 'Color', colors(1,:))

text(x_pos, y_pos(2), text_legend{2}, 'FontSize', fSet.Fontsize_text, 'Color', colors(2,:))

% set ax param

set(gca,'YDir','reverse')

% plot text with stats

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

switch label_drug{idrug}

case 'Dopamine'

x_pos = get_value_range(tmp_x, 0.65);

y_pos = get_value_range(tmp_y, 0.05);

case 'SCH23390'

x_pos = get_value_range(tmp_x, 0.05);

y_pos = get_value_range(tmp_y, 0);

end

h_text(icrit,idrug) = text(x_pos, y_pos, {betaString, pstring_chi}, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3]);

end

Chi(1) 0.42, p = 0.515
Chi(1) 0.82, p = 0.365
Chi(1) 4.71, p = 0.030
Chi(1) 0.74, p = 0.390
Chi(1) 1.03, p = 0.311
Chi(1) 0.80, p = 0.370
Chi(1) 0.69, p = 0.406
Chi(1) 0.80, p = 0.370

% multiple comparison correction, individual models

[p_fdr, p_masked] = FDR(P, 0.05);

plotj_text_emphasise(h_text, p_masked, 'italic', 2);

plotj_text_emphasise(h_text, p_masked, 'bold', 2);

plot gain, mean to variance relationship, separate for unit-types

[colors, labels] = get_colors('att_drug');

[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');

colors_violin = reshape(colors_violin, [4,3]);

labels_violin = reshape(labels_violin, [4,1]);

ncol = length(label_drug);

nrow = 2; % unit types

[fH, fSet] = plotj_initFig('width',20,'height',15,'Journal',plot_conventions);

clear legendText h

gain_table_full = [];

iplot=0;

xval_gain = 1:500;

for idrug = 1:length(label_drug)

for iunitc = 1:length(unique(unitlist.unit_class))

iplot = iplot+1;

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);

% subplot axis

h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

if iunitc==1

plotj_initAx(fSet, 'axlabel', idrug);

else

plotj_initAx(fSet);

end

hold on

if idrug==1

h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);

end

gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');

variable_names = gv.Properties.VariableNames;

for iatt_drug = 1:length(labels_violin)

mus = [gv.(labels_violin{iatt_drug}).mus];

s2s = [gv.(labels_violin{iatt_drug}).s2s];

mus = mus(:);

s2s = s2s(:);

% get binomial fit

rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);

gv_pop = GainVariability;

gv_pop.rhat = rhat;

m2v = gv_pop.pred_variance(xval_gain);

% plot

h(iatt) = scatter(mus,s2s,fSet.MarkerSize*3,colors_violin(iatt_drug,:),'filled','MarkerFaceAlpha',0.5);

plot(m2v.mean,m2v.variance,'Color',colors_violin(iatt_drug,:),'linew',fSet.LineWidth);

xlim([0.8 max(xval_gain)])

ylim([0.8 max(xval_gain)])

plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)

set(gca,'xscale','log','yscale','log')

end

% axis square

% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)

if idrug==length(label_drug)

xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)

end

if iunitc==1

ylabel({[label_drug{idrug} label_drug_ext{idrug}],'Variance (spikes^{2})'},'FontSize',fSet.Fontsize_text)

end

% violin

% init axes

pos = get(h_ax, 'Position');

h_ax_inset = axes('Position', [pos(1)+0.22 pos(2)+0.03 0.16 0.15]) ; % inset

plotj_initAx(fSet);

hold on

% unit selection

idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});

idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);

num_unit_plot = length(unique(att_table.unit(idx_unit)));

fprintf('\t Found %d units \n', num_unit_plot)

% extract smaller table

v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {'gain_log','drug_on','attention','unit','attention_drug_on'}));

% plot violin

v = violinplot(v_att_table.gain_log, v_att_table.attention_drug_on, 'ShowMean', true);

XTickLabel = get(gca,'XTickLabel');

for iv = 1:length(v)

v(iv).ScatterPlot.SizeData = 10;

idx = strcmpi(XTickLabel(iv), labels_violin);

v(iv).ViolinColor = colors_violin(idx,:);

end

set(gca,'XTickLabel',[])

h_ax_inset.XAxis.Visible = 'off';

% plot lines for each unit

[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));

for iv = 1:length(v)

if ~any(isnan(v_att_table.gain_log))

% works if data is complete

data_x(:,iv) = v(iv).ScatterPlot.XData;

data_y(:,iv) = v(iv).ScatterPlot.YData;

else

% catch if nan values are present

idx = strcmpi(XTickLabel(iv), labels_violin);

idx_trial = v_att_table.attention_drug_on==labels_violin{idx};

data_x(:,iv) = iv;

data_y(:,iv) = v_att_table.gain_log(idx_trial);

end

plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);

% fit model

lme_population = fitlme(v_att_table,['gain_log ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding

disp(lme_population.Coefficients)

coefficients = lme_population.Coefficients;

betas = coefficients.Estimate(2:end);

betas_se = coefficients.SE(2:end);

p_val = coefficients.pValue(2:end);

% p_string = get_significance_strings(p_val, 'rounded', 0);

p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});

p_string{end+1} = ['n = ' num2str(num_unit_plot)];

% plot text

YLIM = get(gca,'ylim');

set(gca, 'ylim', YLIM .* [1 1.25])

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.95);

y_pos = get_value_range(tmp_y, 1.3);

h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right');

ylabel('log(\sigma^{2}_G )', 'FontSize', fSet.Fontsize_text)

% [p_fdr, p_masked] = FDR(p_val, 0.05);

[p_masked] = p_val<0.05;

p_masked = [p_masked ; false];

plotj_text_emphasise(h_text, 1, 'italic', p_masked);

plotj_text_emphasise(h_text, 1, 'bold', p_masked);

end

	 Found 14 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate    SE          tStat      DF    pValue        Lower        Upper    
    {'(Intercept)'                        }         -1.6769     0.24109    -6.9553    47    9.5535e-09      -2.1619      -1.1919
    {'drug_on_Drug on'                    }         0.23563    0.074521     3.1619    47      0.002744     0.085713      0.38555
    {'attention_Attend RF'                }        -0.16931    0.074521     -2.272    47      0.027706     -0.31923    -0.019394
    {'drug_on_Drug on:attention_Attend RF'}        0.097133    0.075694     1.2832    47       0.20571    -0.055145      0.24941

	 Found 17 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate     SE          tStat       DF    pValue        Lower       Upper   
    {'(Intercept)'                        }          -1.9154      0.1718     -11.149    63    1.5163e-16     -2.2587      -1.572
    {'drug_on_Drug on'                    }          0.16573    0.045026      3.6809    63    0.00048392    0.075757     0.25571
    {'attention_Attend RF'                }        -0.030159    0.045026    -0.66981    63       0.50543    -0.12014    0.059818
    {'drug_on_Drug on:attention_Attend RF'}        -0.031138    0.045026    -0.69156    63       0.49176    -0.12111    0.058839

	 Found 3 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate     SE          tStat      DF    pValue       Lower         Upper    
    {'(Intercept)'                        }          -1.1245     0.31895    -3.5257    8     0.0077812         -1.86     -0.38901
    {'drug_on_Drug on'                    }        0.0032421    0.029576    0.10962    8       0.91541     -0.064961     0.071446
    {'attention_Attend RF'                }          -0.1411    0.029576    -4.7706    8     0.0014075       -0.2093    -0.072895
    {'drug_on_Drug on:attention_Attend RF'}          0.06162    0.029576     2.0834    8      0.070737    -0.0065831      0.12982

	 Found 11 units 

    Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05

    Name                                           Estimate    SE          tStat      DF    pValue        Lower        Upper   
    {'(Intercept)'                        }        -1.6587      0.34668    -4.7845    37    2.7399e-05      -2.3611    -0.95625
    {'drug_on_Drug on'                    }        0.31791     0.088241     3.6027    37    0.00092093      0.13911      0.4967
    {'attention_Attend RF'                }        -0.1319     0.088241    -1.4948    37       0.14346     -0.31069    0.046895
    {'drug_on_Drug on:attention_Attend RF'}        0.17933     0.087884     2.0405    37      0.048482    0.0012596      0.3574

plot behaviour

mfactor = 1000; %multiplication_factor

colors = get_colors('att_drug');

ncol = 1;

nrow = length(label_drug);

ylim2plot = [250 550];

[fH, fSet] = plotj_initFig('width', 10, 'height', 12, 'Journal',plot_conventions);

clear legendText h

for idrug = 1:length(label_drug)

idx_rec = strcmpi(recordinglist.Drug, label_drug(idrug));

subtightplot(nrow, ncol, idrug, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

plotj_initAx(fSet, 'axlabel', idrug + 1, 'axlabelDisplacement', [0.07, 0.01]);

hold on

title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

% plot errorbar

iattDrug = 0;

[RT_table,RT_mat,colors_mat] = deal([]);

for iatt = 1:2

for idrug_in = 1:2

iattDrug = iattDrug+1;

y = squeeze(mean(RT.RT(idx_rec,idrug_in,abs(idx_attention_cond-2)==iatt),3)) * mfactor;

text_legend{iattDrug} = sprintf('%s, %s', label_attention{iatt}, label_drug_onoff{idrug_in});

n = length(find(idx_rec));

rec = (1:n)';

cond = repmat(iatt,n,1);

drug = repmat(idrug_in,n,1);

RT_table = [RT_table ; table(rec,y,cond,drug,'VariableNames',{'recording','RT','attention','drug'})];

RT_mat = [RT_mat, y];

colors_mat = [colors_mat ; squeeze(colors(iatt,idrug_in,:))'];

end

RT_table.attention = categorical(RT_table.attention);

RT_table.drug = categorical(RT_table.drug);

v = violinplot(RT_mat, text_legend, 'ShowMean', true);

for iv = 1:length(v)

v(iv).ViolinColor = colors_mat(iv,:);

end

% do stats

% check with lme

lme = fitlme(RT_table,['RT ~ 1 + (1|recording)'], 'DummyVarCoding', 'effects');% fit constant

lme2 = fitlme(RT_table,['RT ~ attention + (1|recording)'], 'DummyVarCoding', 'effects'); % fit linear model with attention condition

lme3 = fitlme(RT_table,['RT ~ attention + drug + (1|recording)'], 'DummyVarCoding', 'effects'); % fit linear model with states

lme4 = fitlme(RT_table,['RT ~ attention + drug + attention:drug + (1|recording)'], 'DummyVarCoding', 'effects'); % fit interaction

[BETA,BETANAMES,STATS] = fixedEffects(lme4);

% compare the model fits

comp1 = compare(lme, lme2, 'CheckNesting',true); % compare linear to constant

comp2 = compare(lme2, lme3, 'CheckNesting',true); %

comp3 = compare(lme3, lme4, 'CheckNesting',true); %

[p_fdr, p_masked] = FDR( [comp1.pValue(2) comp2.pValue(2) comp3.pValue(2)], 0.05);

pstring = get_significance_strings([STATS.pValue(2:end)], 'rounded', 0, 'factorstring', {'attention', 'drug', 'interaction'});

pstring{4} = ['n = ' num2str(length(find(idx_rec)))];

ylim(ylim2plot)

% plot text indicating number of units

set(gca, 'Units', 'normalized');

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.1);

y_pos = get_value_range(tmp_y, 0.85);

text(x_pos, y_pos, pstring, 'FontSize', fSet.Fontsize_text)

if idrug==1

ylabel('Reaction time (ms)', 'FontSize', fSet.Fontsize_text)

end

if idrug == length(label_drug)

set(gca, 'xticklabelRotation', 25)

else

set(gca, 'xticklabel', '')

end

pupil timecourse

events = {'stim', 'cue', 'dim'};

xlabel_event = {'Time from stimulus onset (ms)', 'Time from cue onset (ms)', 'Time from first-dimming (ms)'};

xlim2use = {[-200 1000], [-200 1000], [-1000 200]};

ylim2use = [-2.7 0.5];

ncol = length(events); % events

nrow = 2;%length(label_drug)-1;

timewinsize = 200;

timestep = 10;

[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap ./ [1 2];

clear ax h legendText

iplot = 0;

for idrug = 1:nrow

for ievent = 1:length(events)

iplot = iplot+1;

[pupil_timeseries, reclist] = get_population_data(recordinglist, 'pupil_timeseries', path_data, events{ievent});

if nrow==length(label_drug)

idx_rec = strcmpi(reclist.Drug, label_drug(idrug));

else

idx_rec = true(height(reclist),1);

end

nrecs = length(find(idx_rec));

pupil2plot = squeeze(nanmean(pupil_timeseries.pupil(idx_rec,:,:,:),3));

time = pupil_timeseries.time*1000;

timecenters = time(1):timestep:time(end);

p_drug = NaN(1,length(timecenters));

clear legendText

h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);

if ievent==1

plotj_initAx(fSet, 'axlabel', idrug, 'axlabelDisplacement', [0.07, 0.01]);

else

plotj_initAx(fSet);

end

hold on

if (ievent~=1)

axh = get(gca);

axh.YTick = [];

axh.YAxis.Visible = 'off'; % remove y-axis

end

ilegend = 0;

for idrug_in = 1:2

ilegend = ilegend+1;

y = squeeze(nanmean(pupil2plot(:,idrug_in,:),1));

y_err = squeeze(nanstd(pupil2plot(:,idrug_in,:),[],1))/sqrt(size(pupil2plot,1));

h(ilegend) = boundedline(...

time, ...

y, ...

y_err, ...

'alpha', 'cmap', fSet.colors(idrug_in,:));

set(h(ilegend), 'LineWidth', fSet.LineWidth)

text_legend{ilegend} = sprintf('Drug %s', label_drug_onoff{idrug_in});

end

for itc = 1:length(timecenters)

timewin = timecenters(itc) + [-timewinsize/2 timewinsize/2];

idx_time = dsearchn(time(:), timewin(:));

tmp_data = squeeze(nanmean(pupil2plot(:,:,idx_time(1):idx_time(2)),3));

p_drug(itc) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');

end

[p_fdr, p_masked] = FDR( p_drug, 0.05);

plotj_brokenVector(timecenters, ylim2use(1)*0.95, p_masked, 'linewidth', 5, 'color', [0.5 0.5 0.5]);

if iplot==1

hleg = legend(h, text_legend, 'autoupdate', 'off', 'FontSize', fSet.Fontsize_text, 'Box', 'Off');

end

xlim(xlim2use{ievent})

ylim(ylim2use)

plot([0 0], ylim, 'k','linew',1)

if ievent==1

ylabel('Pupil diameter (z-score)', 'FontSize', fSet.Fontsize_text)

end

if idrug==nrow

xlabel(xlabel_event{ievent}, 'FontSize', fSet.Fontsize_text)

end

if ievent~=1

set(gca,'yticklabel',[])

end

if nrow>1

if ievent==2

title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)

end

if ievent==1

text(h_ax, 500, -0.5, sprintf('n=%d',nrecs), 'FontSize', fSet.Fontsize_text)

end

Warning: Ignoring extra legend entries.

Pupil drug/no drug, scatter

colors = get_colors('drug');

events = {'baseline', 'stim', 'cue', 'dim'};

event_title = {'Baseline','Stimulus evoked','Cue evoked','First-dimming'};

ncol = length(events); % events

nrow = 1;

xtext = [22 22];

ytext = [-0.45 -0.32];

[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);

fSet.subplotGap = fSet.subplotGap/1.5;

clear ax h legendText

[P, h_text] = deal(NaN(length(events), length(label_drug)));

for ievent = 1:length(events)

[pupil, reclist] = get_population_data(recordinglist, 'pupil_windows', path_data, events{ievent});

axlim = [];

subtightplot(nrow, ncol, ievent, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)

plotj_initAx(fSet, 'axlabel', ievent+2, 'axlabelDisplacement', [0.05, 0.13]);

hold on

title(event_title{ievent}, 'FontSize', fSet.Fontsize_title)

% significance testing

text_p = cell(length(label_drug),1);

for idrug = 1:length(label_drug)

idx_rec = strcmpi(reclist.Drug, label_drug(idrug));

pupil2plot = squeeze(nanmean(pupil.pupil(idx_rec,:,:),3));

P(ievent,idrug) = compare_means(pupil2plot(:,1), pupil2plot(:,2), 1, 'rank');

text_p{idrug} = get_significance_strings(P(ievent,idrug), 'rounded', 0);

text_legend{idrug} = sprintf('%s (n=%d)', label_drug{idrug}, length(find(idx_rec)));

% mean difference

delta_data = pupil2plot(:,2)-pupil2plot(:,1);

% effect size

d = computeCohen_d(pupil2plot(:,2), pupil2plot(:,1), 'paired');

% print result

fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...

'PUPIL', ...

label_drug{idrug}, ...

'pupil', ...

mean(delta_data), ...

std(delta_data)/sqrt(length(find(idx_rec))), ...

text_p{idrug}, ...

end

pupil2plot = squeeze(nanmean(pupil.pupil,3));

h = plotj_scatter(pupil2plot, ...

'dataIndex', idx_rec+1, ...

'Markersize', markersize, 'markerStyle', {'o','s'}, ...

'MarkerEdgeColor', colors, 'MarkerFaceColor', colors, ...

'MarkerEdgeAlpha', [0.5 0.5] ,'MarkerFaceAlpha', [0.5 0.5]);

% set(gca,'ydir','reverse')

axis square

tmp_x = get(gca,'xlim');

tmp_y = get(gca,'ylim');

x_pos = get_value_range(tmp_x, 0.6);

for idrug = 1:length(label_drug)

y_pos = get_value_range(tmp_y, 0.2-0.1*(idrug-1));

h_text(ievent,idrug) = text(x_pos, y_pos, text_p{idrug}, 'FontSize', fSet.Fontsize_text, 'color', colors(idrug,:));

end

if ievent==1

ylabel({'Drug (z-score)'}, 'FontSize', fSet.Fontsize_text)

end

xlabel({'No drug (z-score)'}, 'FontSize', fSet.Fontsize_text)

if ievent==4

hleg = legend(h, text_legend, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'Location', 'NorthEast');

hleg.Position = hleg.Position + [0.05 0.3 0 0];

end

PUPIL: Dopamine, delta-pupil 0.02 +- 0.04, p = 0.821, Cohens's d = 0.13
PUPIL: SCH23390, delta-pupil 0.01 +- 0.05, p = 0.890, Cohens's d = 0.04
PUPIL: Dopamine, delta-pupil 0.10 +- 0.02, p < 0.001, Cohens's d = 1.09
PUPIL: SCH23390, delta-pupil 0.10 +- 0.03, p = 0.004, Cohens's d = 0.79
PUPIL: Dopamine, delta-pupil 0.08 +- 0.03, p = 0.051, Cohens's d = 0.66
PUPIL: SCH23390, delta-pupil 0.07 +- 0.03, p = 0.031, Cohens's d = 0.62
PUPIL: Dopamine, delta-pupil 0.05 +- 0.03, p = 0.193, Cohens's d = 0.37
PUPIL: SCH23390, delta-pupil 0.06 +- 0.03, p = 0.023, Cohens's d = 0.54

Warning: Ignoring extra legend entries.

[p_fdr, p_masked] = FDR(P, 0.05);

plotj_text_emphasise(h_text, p_masked, 'italic');

plotj_text_emphasise(h_text, p_masked, 'bold');